# Lab 2: Oscillations in Delay Differential Equations # Name: Nirel Gidanian # I worked on this code with: Michelle Steinberg, Talia Shimoni # Please do all of your work for this week's lab in this worksheet. If # you wish to create other worksheets for scratch work, you can, but # this is the one that will be graded. You do not need to do anything # to turn in your lab. It will be collected by your TA at the beginning # of (or right before) next week’s lab. # Be sure to put each exercise in a new cell with an exercise number at the top. # Use enough comments that you and the grader can understand your code. # Label axes on all graphs.

**Function definition:** dde_solve(dde, statevar, delayedvars, history, tmax, timestep)

**Description:** Numerically integrate a delay differential equation

A single-variable delay differential equation is an equation of the form *X*'(*t*) = *f*(*X*, *X*_{τ1}, *X*_{τ2}, ... , *X*_{τn}) where on the right hand side *X* denotes *X*(*t*) and each *X*_{τ} denotes *X*(*t* - τ) for some constant delay τ.

**Inputs:**

dde: The right hand side of the delay differential equation

statevar: The state variable to use (*X* in the example above)

delayedvars: A dictionary whose keys are the symbolic variables that appear in the DDE expression, which represent the value of the state variable some time ago, and whose values are the actual delays. For example, if the DDE is *X*'(*t*) = *X*(*t*) - X(*t* - 5) + *X*(*t* - 17.4) the dde argument might be X - Xtau + Xtau2, in which case this argument would be {Xtau: 5, Xtau2: 17.4}. Note that this means before calling this function, all three of X, Xtau, and Xtau2 should have been declared as symbolic variables: var("X, Xtau, Xtau2")

history: A Sage expression for the history, as a function of t (time)

tmax: Integrate from *t* = 0 to this time

timestep: The size of each time step (the smaller, the more accurate)

**Output:** A list of the values of the state variable (*X* in the example above) at each time step from *t* = 0 to tmax

**Examples:**

A version of Mackey-Glass, using a constant 0 history:

var("X, X_t")

sol = dde_solve(1 - X * X_t^5/(1 + X_t^5), X, {X_t: 8}, history=0, tmax=400, timestep=0.1)

list_plot(zip(srange(0, 400, 0.1), sol), plotjoined=True)

The same, but with an oscillating history function, giving weird transients that give way to a periodic steady state:

var("t, X, X_tau")

eq = 1 - X * X_tau^5/(1 + X_tau^5)

sol = dde_solve(eq, X, {X_tau: 8}, 5*cos(t), 400, 0.1)

list_plot(zip(srange(0, 400, 0.1), sol), plotjoined=True)

%auto def dde_solve(dde, statevar, delayedvars, history, tmax, timestep): # Check validity of delays. if min(delayedvars.values()) < 0: raise ValueError("This function will not work with negative delays. " "Consider consulting a fortune teller instead.") if any(val<timestep for val in delayedvars.values()): raise ValueError("Time step should not be larger than delay.") # Set up variables and delays. delayedvars = delayedvars.items() dde = dde.subs({v: statevar for v, delay in delayedvars if delay == 0}) delayedvars = [(v, delay) for v, delay in delayedvars if delay != 0] allvars = [str(statevar)] + [str(v) for v, delay in delayedvars] delays = [delay for v, delay in delayedvars] # Set up fast functions. dde_func = fast_float(dde, *allvars) history_func = fast_float(history, "t") # Adjust the timestep if necessary mindelay = min(delays) if delays else timestep timestepcorrectionfactor = ceil(timestep / mindelay) timestep /= timestepcorrectionfactor # A function to perform history lookups. def lookup(t): """Does a history lookup at each delay from t, stores result in allvars[1:]""" for i, delay in enumerate(delays): if t - delay <= 0: allvars[i+1] = history_func(t - delay) else: r = (t - delay) / timestep n = floor(r) r -= n allvars[i+1] = result[n]*(1 - r) + result[n + 1]*r # Set up for the first iteration. result = [history_func(0)] lookup(0) for t in sxrange(0, tmax - timestep, timestep): # Compute k1. Note history lookup has already been done. allvars[0] = result[-1] k1 = dde_func(*allvars) # Compute k2. lookup(t + timestep/2) allvars[0] += timestep/2 * k1 k2 = dde_func(*allvars) # Compute k3. Note history lookup has already been done. allvars[0] = result[-1] + timestep/2 * k2 k3 = dde_func(*allvars) # Compute k4. lookup(t + timestep) allvars[0] = result[-1] + timestep * k3 k4 = dde_func(*allvars) # Finally, compute the RK4 weighted average. result.append(result[-1] + (k1 + 2*k2 + 2*k3 + k4)/6 * timestep) return result[::timestepcorrectionfactor]

#EXERCISE 1 var("X") @interact def sigmoid(a=(1,10),b=(1,10),c=(1,10),d=(1,10)): #defines the function and the range of input values f = (a+(b*X^d))/(c^d+X^d) #sets f as the function being plotted g = plot(f, (x,0,10), axes_labels = ["x", "f"]) #sets g equal to the plot of f show(g) #shows the plot

X

Interact: please open in CoCalc

#EXERCISE 2 #increasing "a" will shift the sigmoid function upward along the y-axis #increasing "b" will shift the sigmoid function downward along the y-ais #increasing "c" will shift the sigmoid function downward along the y-axis #increasing "d" affects the sensitivity of the sigmoid and determines whether it is an increasing or decreasing sigmoid

#EXERCISE 3 n = 9 #assigns a value to n k1 = 0.2 #assigns a value to k1 k2 = 0.2 #assigns a value to k2 k3 = 0.2 #assigns a value to k3 var("H,P,G") #defines variables hprime = (1/(1+G^n))-(k1*H) #defines the differential equations for h pprime = H-(k2*P) #defines the differential equations for p gprime = P-(k3*G) #defines the differential equations for g t = srange(0,100,0.1) #creates a list of values from 0 to 100 in step size of 0.1 sol = desolve_odeint([hprime, pprime, gprime], ics=[1,1,1], times = t, dvars = [H,P,G]) #simulates the model list_plot(sol, plotjoined = True, axes_labels = [H,P,G]) #plots the trajectory

(H, P, G)

3D rendering not yet implemented

#runs the same thing as previous cell with different initial values n = 9 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) pprime = H-(k2*P) gprime = P-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, pprime, gprime], ics=[5,2,3], times = t, dvars = [H,P,G]) list_plot(sol, plotjoined = True, axes_labels = [H,P,G]) #plots the trajectory

(H, P, G)

3D rendering not yet implemented

#runs the same thing as previous cell with different initial values n = 9 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) pprime = H-(k2*P) gprime = P-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, pprime, gprime], ics=[2,1,4], times = t, dvars = [H,P,G]) list_plot(sol, plotjoined = True, axes_labels = ["H","P","G"]) #plots the trajectory

(H, P, G)

3D rendering not yet implemented

#Even though the models start with different initial conditions, they all approach a limit cycle attractor.

#EXERCISE 4 #runs the same thing as previous cell with the original initial values and different n - value n = 3 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) pprime = H-(k2*P) gprime = P-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, pprime, gprime], ics=[1,1,1], times = t, dvars = [H,P,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) + list_plot(zip(t,sol[:,2]), color = "green", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - value n = 6 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) pprime = H-(k2*P) gprime = P-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, pprime, gprime], ics=[1,1,1], times = t, dvars = [H,P,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) + list_plot(zip(t,sol[:,2]), color = "green", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - value n = 9 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) pprime = H-(k2*P) gprime = P-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, pprime, gprime], ics=[1,1,1], times = t, dvars = [H,P,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) + list_plot(zip(t,sol[:,2]), color = "green", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - value n = 12 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) pprime = H-(k2*P) gprime = P-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, pprime, gprime], ics=[1,1,1], times = t, dvars = [H,P,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) + list_plot(zip(t,sol[:,2]), color = "green", plotjoined = True) #plots the time series

(H, P, G)

#for the value n = 12, there are persistent oscillations since the amplitude of the oscillations remains the same

#EXERCISE 5 #runs the same thing as previous cell with the original initial values and different n - values. The difference is that the P variable is no longer accounted for, so there is no time delay and only two differential equations are plotted and simulated. n = 6 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,G") hprime = (1/(1+G^n))-(k1*H) gprime = H-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, gprime], ics=[1,1], times = t, dvars = [H,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - values. The difference is that the P variable is no longer accounted for, so there is no time delay and only two differential equations are plotted and simulated. n = 9 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) gprime = H-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, gprime], ics=[1,1], times = t, dvars = [H,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - values. The difference is that the P variable is no longer accounted for, so there is no time delay and only two differential equations are plotted and simulated. n = 12 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) gprime = H-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, gprime], ics=[1,1], times = t, dvars = [H,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - values. The difference is that the P variable is no longer accounted for, so there is no time delay and only two differential equations are plotted and simulated. n = 15 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) gprime = H-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, gprime], ics=[1,1], times = t, dvars = [H,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) #plots the time series

(H, P, G)

#runs the same thing as previous cell with the original initial values and different n - values. The difference is that the P variable is no longer accounted for, so there is no time delay and only two differential equations are plotted and simulated. n = 100 k1 = 0.2 k2 = 0.2 k3 = 0.2 var("H,P,G") hprime = (1/(1+G^n))-(k1*H) gprime = H-(k3*G) t = srange(0,100,0.1) sol = desolve_odeint([hprime, gprime], ics=[1,1], times = t, dvars = [H,G]) list_plot(zip(t,sol[:,0]), color = "blue", axes_labels = ["time", "state variables"], plotjoined = True) + list_plot(zip(t,sol[:,1]), color = "red", plotjoined = True) #plots the time series

(H, P, G)

#None of the oscillations, no matter what value of n, have persistent oscillations because the two criteria for oscillations, time delay, and negative feedback are not present, like they are in the model simulated in exercise 4. In this case, there is only a negative feedback loop, but no time delay because the middle-man, the pituitary gland, is eliminated from the model.

#EXERCISE 6

#EXERCISE 7 menu = {"ice cream":4.95, "cookies":2, "candy":1} #creates a dictionary

menu["cookies"] #tests the dictionary

2

#EXERCISE 8 var("X,Y") #defines the varibales L = 6 #assigns a value to X Vmax = 80 #assigns a value to Vmax a = 0.2 #assigns a value to a n = 5 #assigns a value to n h = 400 #assigns a value to h tau = 0.2 #assigns a value to tau prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) #calls the entire equation prime t = srange(0, 100, 0.01) #creates a list starting from 0 and going to 100 in step size of 0.1 sol = dde_solve(prime,X,{Y:0.2},0,100,0.01) #simulates the model list_plot(zip(t,sol), axes_labels = ["Time", "X"], plotjoined = True) #plots the time series

(X, Y)

#EXERCISE 9 var("X,Y") #identifies variables graph = [] #creates an empty list t = srange(0, 100, 0.01) #creates a list starting at 0 to 100 in step size 0.1 for i in sol: #creates a loop w index i going through values of the simulation prime = (Vmax*(i^n))/(h+i^n) #assigns the equation to prime graph.append(prime) #adds the values of prime to the end of graph

(X, Y)

#EXERCISE 10 list_plot(zip(t,graph), axes_labels = ["Time", "Breathing Rate"], plotjoined = True) #plots the trajectory

#the graph eventually settles at a steady oscillation with an amplitude much larger than the amplitude of the graph modelled in exercise 8. Both these graphs represent a healthy breathing rate.

#EXERCISE 11 #does the same thing as exercise 8 with a different value for tau var("X,Y") L = 6 Vmax = 80 a = 0.2 n = 5 h = 400 tau = 0.5 prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) t = srange(0, 100, 0.01) sol = dde_solve(prime,X,{Y:tau},0,100,0.01) list_plot(zip(t,sol), color = "blue", axes_labels = ["Time", "X"], plotjoined = True)

(X, Y)

#does the same thing as exercise 8 with a different value for tau var("X,Y") L = 6 Vmax = 80 a = 0.2 n = 5 h = 400 tau = 1 prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) t = srange(0, 100, 0.01) sol = dde_solve(prime,X,{Y:tau},0,100,0.01) list_plot(zip(t,sol), color = "blue", axes_labels = ["Time", "X"], plotjoined = True)

(X, Y)

#does the same thing as exercise 8 with a different value for tau var("X,Y") L = 6 Vmax = 80 a = 0.2 n = 5 h = 400 tau = 1.5 prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) t = srange(0, 100, 0.01) sol = dde_solve(prime,X,{Y:tau},0,100,0.01) list_plot(zip(t,sol), color = "blue", axes_labels = ["Time", "X"], plotjoined = True)

(X, Y)

#All values cause oscillation because so long as there is a delay in the time the brain receives the information, an oscillation will occur. The magnitude of the time delay effects the wavelength of the oscillation, so for some values of tau that are extremely small, the oscillation will be hard to detect. The only limitation on the time delay value, tau, is that it must be larger than the step size.

#EXERCISE 12 #turns the code from exercise 8 into an interative to manipulate tau. @interact def tauchange(tau=(0.02,1.5)): var("X,Y") L = 6 Vmax = 80 a = 0.2 n = 5 h = 400 prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) t = srange(0, 100, 0.01) sol = dde_solve(prime,X,{Y:tau},0,100,0.01) k = list_plot(zip(t,sol), color = "blue", axes_labels = ["Time", "X"], plotjoined = True) show(k)

Interact: please open in CoCalc

#EXERCISE 13 #0.17688 is the approximate minimum value of tau necessary to create a persistent oscillation.

#EXERCISE 14 #the parameter n controls steepness.

#EXERCISE 15 #turns the code from exercise 12 into an interative to manipulate tau and n. @interact def tauchange(tau=(0.02,1.5), n=(0.1,10)): var("X,Y") L = 6 Vmax = 80 a = 0.2 h = 400 prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) t = srange(0, 100, 0.01) sol = dde_solve(prime,X,{Y:tau},0,100,0.01) k = list_plot(zip(t,sol), color = "blue", axes_labels = ["Time", "X"], plotjoined = True) show(k)

Interact: please open in CoCalc

#EXERCISE 16 #same code as exercise 15 @interact def tauchange(tau=(0.02,1.5), n=(0.1,10)): var("X,Y") L = 6 Vmax = 80 a = 0.2 h = 400 prime = L-(((Vmax*(Y^n))/(h+Y^n))*a*X) t = srange(0, 100, 0.01) sol = dde_solve(prime,X,{Y:tau},0,100,0.01) k = list_plot(zip(t,sol), color = "blue", axes_labels = ["Time", "X"], plotjoined = True) show(k)

Interact: please open in CoCalc

#the minumum n value that results in persistent oscillation is approximatley 4.9906

#EXERCISE 17 tlist = [0.50544,0.5676,0.759,1.00864,1.355] nlist = [2.9314,2.7928,2.476,2.2186,2.0008]

#EXERCISE 18 list_plot(zip(tlist,nlist), axes_labels = ["tau","sensitivity"], plotjoined = True) #plots tau and sensitivity against each other in pairs.

#as tau increases, the minimum sensitivity for persistent oscillation decreases, which shows that they are inversly related.