I would like to solve a system of network-based differential equations using python 3.6

Question

I would like to solve a system of network-based differential equations using python 3.6. The system of equations is as follows:

dx_i/dt = omega_i  - epsilon_i * y_i * x_i - mu_i * x_i,  

dy_i/dt = epsilon_i * y_i * x_i - zeta_i * y_i - rho_i * y_i * z_i,

dv_i/dt = c_i * y_i - gamma_i * v_i + \sum_{i neq j} beta_{ij} * v_i * x_i,

dz_i/dt = k_i * y_i * z_i - delta_i * z_i,

where beta_{ij} = beta (1 - sigma_{ij}) * exp(- alpha|i-j|) 

i = 1,2,3,...,N

I have written the code below in an attempt to solve the system of differential equations in a spatial network.

from jitcode import jitcode, y
import numpy as np
import sympy

#import symengine
import matplotlib.pyplot as plt
#from scipy.integrate import ode
from numpy.random import uniform

n = 10
alpha = 0.05
#beta = 0.1

beta = uniform(0.01,3.0,n)
beta.sort()
mu = uniform(0.01,3.0,n)
mu.sort()
epsilon = uniform(0.01,3.0,n)
epsilon.sort()
pi = uniform(0.01,3.0,n)
pi.sort()
gamma = uniform(0.01,3.0,n)
gamma.sort()
omega = uniform(0.01,3.0,n)
omega.sort()
zeta = uniform(0.01,3.0,n)
zeta.sort()
rho = uniform(0.01,3.0,n)
rho.sort()
k = uniform(0.01,3.0,n)
k.sort()
c = uniform(0.01,3.0,n)
c.sort()


# Knonecker delta

M = np.einsum('ij-> ij',np.eye(n,n))

print(M)

# Adjacency matrix

A = beta * M * sympy.exp(-alpha) 

print(A)


def f():
    for i in range(n):

    coupling_sum = sum(y(i+2) * y(i) for j in range(n) if A[i, j]
    )
    yield omega[i] - epsilon[i] * y(i+2) * y(i) - mu[i] * y(i)
    yield epsilon[i] * y(i+2) * y(i) - zeta[i] * y(i+1) - rho[i] * y(i+1)* y(i+3)
    yield c[i] * y(i+1) - gamma[i] * y(i+2) + coupling_sum
    yield k[i]* y(i+1) * y(i+3) - delta[i] *y(i+3)

#integrate
#---------------

initial_state = np.random.random(n)

ODE = jitcode(f,n=n)
ODE.set_integrator("dopri5", atol=1e-6,rtol=0)
initial = np.linspace(0.1,0.4,num=n)

ODE.set_initial_value(initial_state,time=0.0)
#data structure: x[0], w[0], v[0], z[0], ..., x[n], w[n], v[n], z[n]
data = []
data = np.vstack(ODE.integrate(T) for T in range(0, 100, 0.1))
print(data)
fig = plt.figure()

I am expecting to get solutions for the four differential equations and some simulations to represent the equations. The error message that Iam getting says "RuntimeError: Not Implemented"

does it comes from jitcode?which line is raising the error? also you seem to have indentation error in f — Benoit de Menthière, Aug 11 '19 at 04:49
Did you mean the linked equation to be `dv_i/dt = c_i * y_i - gamma_i * v_i - \sum_{j neq i} beta_{ij} * v_j * x_j`? The implementation of this equation can not really be connected to what the equation says, neither can be used to make sense of the other. — Lutz Lehmann, Aug 11 '19 at 08:56
The error message is coming from the line before print(data) code. i.e. data = np.vstack.... The message now says TypeError: 'float' object cannot be interpreted as an integer. The coupling_sum is the invading viral particles for a population of n individuals that are connected in a network, thus it is added to the viral population on dv_i / dt — BLESSINGS MUFOYA, Aug 12 '19 at 19:10
I used the adjacency matrix A[i j] to model beta_{ij}. How can i do this better? — BLESSINGS MUFOYA, Aug 12 '19 at 19:30
Please comment further on what your system represents. I read from your formulas that in each cell the basis population is `x`, which is constantly created and has a die-off rate. Whenever an `x` and and `y` meet, there is some chance that the `x` gets converted into an `y`. The same with `y` and `z`. On top of that there is the `v` population that has a creation by spontaneous conversion from `y`, dies at some rate and is created when a local `v` meets a cell-foreign `x`. — Lutz Lehmann, Aug 13 '19 at 06:42
Please do not edit your question to a different one after it has been answered. I rolled back your last edit to keep things consistent. If you have a new question on plotting, ask it separately, preferably using some simple array as an example and not something created by a complicated ODE. — Wrzlprmft, Aug 15 '19 at 06:03

score 1 · Answer 1 · answered Aug 13 '19 at 06:19

Your translation of the model into code has index problems. You addressed this once in the sum calculation, but then never again. To help with that, define helper functions

def X(i): return y(4*i)
def Y(i): return y(4*i+1)
def V(i): return y(4*i+2)
def Z(i): return y(4*i+3)

Then you can express the generator of the symbolic right sides as

def f():
    for i in range(n):

        coupling_sum = V(i) * sum(beta[i,j]*X(j) for j in range(n) if j!=i )
        yield omega[i] - epsilon[i] * Y(i) * X(i) - mu[i] * X(i)
        yield epsilon[i] * Y(i) * X(i) - zeta[i] * Y(i) - rho[i] * Y(i)*Z(i)
        yield c[i] * Y(i) - gamma[i] * V(i) + coupling_sum
        yield k[i] * Y(i) * Z(i) - delta[i] * Z(i)

with an appropriate definition of beta[i,j]. The code you wrote is strange and does not match the formula. In the formula, beta is first a constant and then a matrix. In the code, beta is an array. This is quite incompatible.

In the call to compile the function, you also should give the correct dimension of the state space, you have n components of x,y,v,z, making 4*n components in total.

initial_state = np.random.random(4*n)
ODE = jitcode(f,n=4*n)

With the correct dimension of the state space arguments, the integration calls probably will go through.

I would like to solve a system of network-based differential equations using python 3.6

1 Answers1