Delay differential equations in python

Question

I have this control problem, where my goal is to select values of D such that Hb(t) stays within a certain range, see the equations here (if it helps I'm trying to reproduce the simulation discussed in this paper).

I'm not confident in my attempt to solve this problem- specifically the terms where past state values of E and P are needed. My effort was to calculate them manually and pass as arguments to my integrating function. But this approach:

• has repeated code and is difficult to understand
• is incorrect I think because as after solving, I get Hb(t) values that I would not expect even after changing many parameters values.

Is there a better way to solve these differential equations with time related states?

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def erythropoiesis_func(s, t, d, term1, term2, term3, term4):
    """
    Inputs
    d     = drug dose (ug/kg per week)
    term1 = avg of past Tp days
    term2 = E_{tot}(t - (Tp + Tm))
    term3 = P(t - (Tp + Tm))
    term4 = sum of past Tj days times normal
    """

    """
    States of model s(E,P,R):

    E: exogenous plasma concentration of EPO
    P: bone marrow concentration of progenitor cells and precursor cells
    R: concentration of red blood cells (RDCs)

    """
    E = s[0]                
    P = s[1]                
    R = s[2]

    # Parameters to vary by each trajectory (patient) but constant for now
    Ep = .36  #np.random.normal(.3588, .0753, 1)
    Cr = .14 #np.random.normal(.1372, .0520, 1)
    Cp = .21 #np.random.normal(.2014, .0640, 1)

    # Constants:
    Vd  = 52.4            # Volume distribution
    E50 = 100 / Vd        # half-maximal effective concentration (drug potency)
    Tp  = 9               # P cell life in days
    Tm  = 4               # Delay from P compartment to R compartment in days
    Tr  = 70              # R cell life in days


    Etot = (d / Vd) + Ep

    Spmr = np.random.normal(Tr , 30, size=Tr).sum()


    # Compute sdot:
    dsdt = np.empty(len(s))
    dsdt[0] = (24/25)*np.log(2)*Etot                      # E'(t)
    dsdt[1] = Cp*(Etot/(E50 + Etot))*P - Cp*term1         # P'(t)
    dsdt[2] = Cr*(term2/(E50 + term2))*term3 - (Cr/Spmr)*term4

    return dsdt

# Initial Condition for the Control
d0 = .5

# Initial Conditions for the States (E, P, R)
s0 = np.array([(d0/52.4), 1., 1.])

# Constants
Tp  = 9
Tm  = 4               
Tr  = 70
E50 = 100 / 52.4
MCH = 2.4

# Time Interval (days)
tt = np.array(range(Tp+Tm+Tr))

# Pre fill vectors
E = np.ones(len(tt)) * s0[0]
P = np.ones(len(tt)) * s0[1]
R = np.ones(len(tt)) * s0[2]
H = MCH * R

d = np.random.uniform(low=0.0, high=1.0, size=len(tt))

# Simulate
term1_vals = []
term4_vals = []
for t in range(len(tt)-1):
    # time to feed to solver
    ts = [tt[t],tt[t+1]]

    # generating term1
    temp1 = []
    for tj in range(Tp):
        if tj <= t:
            temp1.append( (E[t - tj] / (E50 + E[t - tj] )) * P[t - tj] )
    term1_vals.append(np.mean(temp1))

    # collecting term2 and term3 if they exist
    if (Tp + Tm) <= t:
        term2 = E[t-(Tp + Tm)]
        term3 = P[t-(Tp + Tm)]
    else:
        term2 = 0
        term3 = 0

    # generating term4
    temp4 = []
    for tj in range((Tp+Tm),(Tp+Tm+Tr)):
        if tj <= t:
            temp4.append( np.random.normal(Tr , 30, size=1) * (E[t - tj] / (E50 + E[t - tj] )) * P[t - tj] )
    term4_vals.append(np.sum(temp4))

    # calling solver
    s = odeint(erythropoiesis_func,s0,ts,args=(d[t+1],term1_vals[t],term2,term3,term4_vals[t]))

    #storing values
    E[t+1] = s[-1][0]
    P[t+1] = s[-1][1]
    R[t+1] = s[-1][2]
    H[t+1] = MCH * s[-1][2]
    s0 = s[-1]