Reputation: 21
Solving Coupled ODE's in Python - too much/little oscillation
So I've been trying to model a system of coupled, first order differential equations that I found in a paper. They are of biological nature, so negative solutions for example don't make any sense.
Since I know what the final figure should look like from the paper (see below), I know that my solution is plain wrong.
Figure from paper:
In particular, there are two scenarios I have managed to create so far with my code:
- Default solver: Heavy oscillation of no and np, while nc and ns stay constant - should not be the case.
- BDF Solver: Immediate "steady-state" (not the actual steady state of the system though) after 0 seconds. Basically useless.
This is my code (scenario 1), I have decided to only use three out of the four DE's to make room for the constraint that ntitin should stay constant.
import numpy as nump
# Constants
kp = 0.07 # M^-1 s^-1
kr = 6 # s^-1
ks = (10**(-8)) - (10**(-6)) # s^-1
kdns = 0.002 # s^-1
ATP = 7.2*10**(-3) # M
w0 = 10**6 # s^-1
kb = 1.3806452*10**(-23) # J/K
T = 310 # K
u0 = 1*10**(-9) # m
umax = 28*10**(-9) # m
delta_G = 35*kb*T
ntitin = 4.7 * 10**(-6) # in molar (M)
# Differential equations
def equations(t, y):
nc, no, np = y
# Calculate ns using the constraint
ns = ntitin - (no + nc + np)
# Calculate Eplus and Eminus for the current force value
f = force_input(t) # Call the previously defined force input function
Eplus = delta_G - f * u0
Eminus = f * (umax - u0)
# Calculate rate constants for the current force value
kplus = w0 * nump.exp((-Eplus) / (kb * T))
kminus = w0 * nump.exp((-Eminus) / (kb * T))
# System of differential equations (without dns_dt since ns given through constraint )
dnc_dt = -kplus * nc + kminus * no
dno_dt = kplus * nc - kminus * no - kp * no * ATP + kr * np + kdns * ns
dnp_dt = kp * no * ATP - kr * np - ks * np
#dns_dt = ks * np - kdns * ns
return [dnc_dt, dno_dt, dnp_dt]
# Initial conditions
no0 = 0.0006*ntitin
np0 = 0.0025*ntitin
ns0 = 10e-12
nc0 = ntitin-no0-np0-ns0 #assuming ns starts at 10e-12 and ntitin needs to stay constant
y0 = [nc0, no0, np0]
# IF NS STARTS AT 10E-6, ALL OF THE SUDDEN no INCREASES DRAMATICALLY OVER THE TIME PERIOD
#, BUT NC, NP AND NS STAY STEADY
# Time span for the solver
t_span = [-1, 3 * 260] # 3 full cycles of the input
t_values = nump.linspace(t_span[0], t_span[1], 1000)
# Solve the system
solution = solve_ivp(
equations, t_span, y0, t_eval = t_values)
# Calculate ns for each time point using the constraint
ns_values = ntitin - (solution.y[0] + solution.y[1] + solution.y[2])
sum_values = (ns_values + solution.y[1] + solution.y[1])/ntitin
# Plotting the results
#plt.plot(solution.t, solution.y[0]/ntitin, label='nc')
plt.plot(solution.t, solution.y[1]/ntitin, label='no')
plt.plot(solution.t, solution.y[2]/ntitin, label='np')
plt.plot(solution.t, sum_values, label='ns+np+no', linestyle='--') # Dashed line for ns
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Fraction of open/phosphorylated Titin')
plt.title('System Dynamics Over Time with Constraint')
plt.show()`
The outcome is this:
The force input is correct - I have verified that with additional plots.
Edit: Here is the force input function:
def force_input(t):
"""
Generates a force input as a function of time with 10 on-off repetitions followed by a rest.
:param t: Time in seconds
:return: Force in N at the given time
"""
high_force = 20 * 10**(-12) # N for 'on' state
low_force = 4.5 * 10**(-12) # N for 'off' state and rest
# Period for one 'on-off' cycle (20 seconds), and full period for the entire set (260 seconds)
on_off_period = 20
full_set_period = 260
# Determine the position in the current set
position_in_set = t % full_set_period
if t<0:
return low_force
# Check if we are in the rest period
if position_in_set > 10 * on_off_period:
return low_force
else:
# If not in rest, determine if we are in 'on' or 'off' within the on-off period
position_in_on_off = position_in_set % on_off_period
return high_force if position_in_on_off < 10 else low_force
Also, if I plot the concentrations at the first 10 timepoints, here is the outcome:
nc: [0.99689787 0.99689787 0.99689785 0.99689779 0.99689774 0.99689769
0.99689764 0.99689758 0.99689753 0.99689748]
no: [ 0.0006 0.19058381 -0.33508352 -0.02302532 -0.01392478 -0.14079521
0.01207641 -0.07641563 0.0918557 0.04086486]
np: [ 0.0025 -0.18748378 0.3381835 0.0261254 0.01702492 0.14389538
-0.00897615 0.07951593 -0.08875532 -0.03776443]
ns: [1.00000000e-11 9.83700948e-12 1.02291693e-11 9.97153525e-12
9.94892783e-12 1.00318165e-11 9.89773534e-12 9.95090262e-12
9.80492439e-12 9.82905055e-12]
Seems like no and np just compensate each other through oscillation even into the negative space.
Upvotes: 2
Views: 148
Answers (1)
Reputation: 6745
I'm not sure that I could call this an answer - just a work in progress. However.
Your equations are very stiff. (They contain exponential-decay-like terms.) Basically, you need a solver that can deal with stiff equations; I suggest either 'Radau' or 'BDF'. Usually I would be looking for k.dt << 1, where k is the largest rate constant. To force the timestep dt use the max_step= option. (It needs to be quite low).
You should initialise nc, n0, np and ns with something satisfying the "steady-state" solution (time derivative zero and no forcing). Your current values don't.
I'm not sure that your rate constants are correct. The paper (or its supplementary material) suggest that kp.ATP is roughly 5.kr, so I've taken the liberty of enforcing that.
With the forcing that you use, the rate constants kplus and kminus are just too small compared to the other rate constants to do much (though I'm not entirely sure that the given value of kr is correct in the paper; it looks too large.) I've arbitrarily increased w0 for now, but it's a little unsatisfactory.
I don't think that enforcing the constraint on nc+np+no+ns is particularly helpful, here. It leads to unnatural asymmetry amongst the equations, which should satisfy the constraint perfectly well themselves (just add them up).
At the end of the day, I think your best bet is to email the Cambridge academics who wrote the paper and ask them precisely which parameters they used to create their figure: it is VERY sensititive to the numbers used.
import numpy as nump
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# Constants
kp = 0.07 # M^-1 s^-1
kr = 6 # s^-1
ks = 1.0e-7 # s^-1
kdns = 0.002 # s^-1
ATP = 7.2e-3 # M
w0 = 1.0e6 # s^-1
kb = 1.3806452e-23 # J/K
T = 310 # K
u0 = 1.0e-9 # m
umax = 28.0e-9 # m
delta_G = 35 * kb * T
ntitin = 4.7e-6 # in molar (M)
# OVER-RIDDEN VALUES - please see paper and supplementary material
kp = 5 * kr / ATP # Paper says that ATP.kp is about 5 kr
w0 = 100 * w0 # Otherwise forcing is just too small
def force_input(t):
"""
Generates a force input as a function of time with 10 on-off repetitions followed by a rest.
:param t: Time in seconds
:return: Force in N at the given time
"""
high_force = 20e-12 # N for 'on' state
low_force = 4.5e-12 # N for 'off' state and rest
# Period for one 'on-off' cycle (20 seconds), and full period for the entire set (260 seconds)
on_off_period = 20
full_set_period = 260
# Determine the position in the current set
position_in_set = t % full_set_period
if t<0:
return low_force
# Check if we are in the rest period
if position_in_set > 10 * on_off_period:
return low_force
else:
# If not in rest, determine if we are in 'on' or 'off' within the on-off period
position_in_on_off = position_in_set % on_off_period
return high_force if position_in_on_off < 10 else low_force
# return high_force if int( t / on_off_period ) % 2 == 1 else low_force # Simpler alternative
# Differential equations
def equations(t, y):
nc, no, np, ns = y
# Calculate Eplus and Eminus for the current force value
f = force_input(t) # Call the previously defined force input function
Eplus = delta_G - f * u0
Eminus = f * (umax - u0)
# Calculate rate constants for the current force value
kplus = w0 * nump.exp( -Eplus / (kb * T))
kminus = w0 * nump.exp( -Eminus / (kb * T))
# print( kminus, kplus, kp * ATP, kr, kdns ) # Check the size of rate constants
# System of differential equations
dnc_dt = -kplus * nc + kminus * no
dno_dt = kplus * nc - kminus * no - kp * ATP * no + kr * np + kdns * ns
dnp_dt = kp * ATP * no - kr * np - ks * np
dns_dt = - kdns * ns + ks * np
return [dnc_dt, dno_dt, dnp_dt, dns_dt]
# Initial conditions - SATISFY STEADY STATE WITH NO FORCING
no0 = 0.0006 * ntitin # Fix this value
np0 = no0 * kp * ATP / ( ks + kr ) # Then everything else follows for the steady-state solution
ns0 = np0 * ks / kdns
nc0 = ntitin - no0 - np0 - ns0 # To all intents and purposes nc0 = ntitin
y0 = [ nc0, no0, np0, ns0 ]
# Time span for the solver
t_span = [-50, 3 * 260 ] # steady, followed by 3 full cycles of the input
t_values = nump.linspace(t_span[0], t_span[1], 100000 )
# Solve the system
solution = solve_ivp( equations, t_span, y0, t_eval=t_values, method='Radau', max_step=0.02 )
sum_values = ( solution.y[1] + solution.y[2] + solution.y[3] ) / ntitin
# Plotting the results
plt.plot(solution.t, solution.y[1]/ntitin, label='no')
plt.plot(solution.t, solution.y[2]/ntitin, label='np')
plt.plot(solution.t, sum_values, label='ns+np+no', linestyle='--') # Dashed line for ns
plt.legend()
plt.xlabel('Time (s)')
plt.ylabel('Fraction of open/phosphorylated Titin')
plt.title('System Dynamics Over Time')
plt.show()
Output (note: quite slow):
Upvotes: 1
Related Questions
- Why is reading lines from stdin much slower in C++ than Python?
- Errors while solving ODE's python
- Fixing Delayed Time Error "NDSolveValue::rdelay" in Wolfram Mathematica for a Dengue Epidemic Model with Delay Differential Equations
- Getting derivative value from JiTCDDE instance
- Unexpected result when solving a system of coupled ode's implicitely
- Solving Second Order ODE's in Matlab
- sympy system of differential equations: Not Implemented