Writing differential equations (to be used with solve_ivp) from a CSV file

I wrote a module called JiTCODE which allows you to specify a differential equation in symbolic form. This is convenient for you since it mostly takes care of translating your equations to a function. Also, during the actual integration, your function evaluation will be very efficient as it is hard-coded in C. This probably matters if you have a hundred and fifty equations.

For your specific problem, I would do it like this:

from jitcode import jitcode,y
import numpy as np
from pandas import read_csv

# Dictionary describing the RHS of the ODE:
f = {}
def add_to_equation(dynvar,to_be_added):
    try:
        f[dynvar] += to_be_added
    except KeyError:
        f[dynvar]  = to_be_added

# Dictionary mapping identifiers to JiTCODE’s dynamic variables
dynvars = {}
def get_dynvar(name):
    try:
        return dynvars[name]
    except KeyError:
        new_dynvar = y(len(dynvars))
        dynvars[name] = new_dynvar
        return new_dynvar

# Parsing the file
for source_name,target_name,rate in read_csv("flows.csv").values:
    rate = float(rate)
    source = get_dynvar(source_name)
    target = get_dynvar(target_name)

    add_to_equation(target, rate*source)
    add_to_equation(source,-rate*source)

initial_condition = np.random.random(len(dynvars))

# Your test case:
Liver_1 = get_dynvar("Liver_1")
Blood_5 = get_dynvar("Blood_5")
assert f[Liver_1] == 0.462*Blood_5-0.0009242*Liver_1-0.045286*Liver_1

# Setup integrator
ODE = jitcode(f)
ODE.set_integrator("dopri5")
ODE.set_initial_value(initial_condition,0.0)

# Actual integration
times = np.linspace(0,100,50)
for time in times:
    print(ODE.integrate(time))

Some notes:

• Including different, non-linear components should be straightforward.

• I here print the output, but you store further process it in different ways.

• You can also use ODE.y_dict to access the solution as a dictionary using your identifiers (via dynvars). This way you never need to directly deal with the numbering of dynamical variables.

• Basic simplifications like

  -0.000924*Liver_1 - 0.045286*Liver_1 = -0.04621*Liver_1
  

  will be done on the fly under the hood by SymEngine (the symbolic backbone).

• JiTCODE uses scipy.integrate.ode or solve_ivp under the hood.

I'm not sure how .solve_ivp() works, but if you just want the string representation of the equation I would do the following

import pandas as pd

df = pd.read_csv("test.csv")  # create DataFrame

def get_dt(some_variable):
    equation = ["d{}/dt =".format(some_variable)]  # store str representation of our eqn in a list
    variables = []
    for index, row in df.iterrows():  # iterate over all rows to check for a match
        if (row["From"] == some_variable) or (row["To"] == some_variable):  # if variable is found in row
            if len(variables):
                variables.append("- {}*{}".format(row["Rate"], row["From"]))  # add it to our equation
            else:
                variables.append("{}*{}".format(row["Rate"], row["From"]))  # the first number is positive for some reason
    equation.extend(variables)  # combine left and right hand sides into one list
    return " ".join(equation)  # combine our list or strings and return a nicely formatted single string

eqn = get_dt("Liver_1")
print(eqn)

Output

dLiver_1/dt = 0.462*Blood_5 - 0.0009242*Liver_1 - 0.045286*Liver_1

Additionally if you want to get all the equations in one go

all_vars = set(df["From"])  # make a "list" of the From variables that doesn't contain duplicates
all_vars.update(set(df["To"]))  # add in all To variables still without duplicates
for var in all_vars:
    eqn = get_dt(var)
    print(eqn)

Which gives us the output

dT-bone-S/dt = 0.12474*Blood_5
dT-bone-V/dt = 0.013859999999999999*Blood_5
dUB-cont/dt = 0.0154*Blood_5 - 3.5*Blood_4 - 0.017329*Kidneys_1
dRC-cont/dt = 0.01155*Blood_5
dSI-cont/dt = 0.0009242*Liver_1
dOther_3/dt = 300.0*Blood - 0.099*Other_3 - 28.95*Blood_4
dLiver_2/dt = 0.045286*Liver_1
dOther_5/dt = 0.0231*Blood_5 - 0.0001266*Other_5
dLiver_1/dt = 0.462*Blood_5 - 0.0009242*Liver_1 - 0.045286*Liver_1
dOvaries/dt = 8.47e-05*Blood_5
dKidneys_1/dt = 0.0077*Blood_5 - 0.017329*Kidneys_1
dKidneys_2/dt = 0.000385*Blood_5 - 0.0001266*Kidneys_2
dBlood/dt = 300.0*Blood - 700.0*Blood
dOther_4/dt = 0.018511*Blood_5 - 0.0013859999999999999*Other_4
dBlood_5/dt = 700.0*Blood - 0.462*Blood_5 - 0.08778*Blood_5 - 0.00462*Blood_5 - 0.12474*Blood_5 - 0.013859999999999999*Blood_5 - 0.0154*Blood_5 - 0.0077*Blood_5 - 0.000385*Blood_5 - 0.01155*Blood_5 - 0.0002695*Blood_5 - 8.47e-05*Blood_5 - 0.018511*Blood_5 - 0.0231*Blood_5 - 0.099*Other_3 - 67.55*Blood_4
dBlood_4/dt = 3.5*Blood_4 - 67.55*Blood_4 - 28.95*Blood_4 - 0.0001266*Kidneys_2 - 0.0013859999999999999*Other_4 - 0.0001266*Other_5
dTestes/dt = 0.0002695*Blood_5
dC-bone-S/dt = 0.08778*Blood_5
dC-bone-V/dt = 0.00462*Blood_5

I think that's enough to get you going so you can figure it out, but it's probably not exactly what you're looking for because of odd lines like:

dBlood/dt = 300.0*Blood - 700.0*Blood

Which doesn't seem to make much sense to me.