Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import numpy as numpy
- from scipy import integrate
- from matplotlib.pylab import *
- def calcium(t, x):
- """
- returns
- """
- n = 1
- #PTH constants
- SetPoint_PTH_50 = 46.8 #mg/L
- PTHmax = 127*10e-6 #mg/L
- dpth = log(2)/(4.0/(60.0*24.0)) # /days
- pth_mean = 35.8e-6 # mg/L
- #CTN constants
- SetPoint_CT = 53.2 #mg/L
- dctn = log(2)/(10.0/(24.0*60.0)) # /days
- CTNmax = 1.0e-5 # mg/L
- #DHCC
- DHCCmean = 2.163e-2 #mg/L
- ddhhc = log(2)/(5.0/24.0) #/days
- #calcium equation constants
- K_PTH = PTHmax*dpth #mg/L*day
- K_CTN = CTNmax*dctn #mg/L*day
- K_DHCC = (DHCCmean*ddhhc)/pth_mean #/day
- K1 = 10000 #PTH
- K2 = 1 #CTN
- K3 = (9.24e-5)*K2 - 0.00165*K1 + 2323.62 #dthh
- K4 = 10.0 #mg/L * day
- K5 = 0.1 #mg/L*day
- #Assign some variables for convenience of notation
- ca = x[0]
- pth = x[1] #mean PTH 35.8 pg/mL
- dthh = x[2] #mean DTHH 51.5 nmol/L
- ctn = x[3] #mean 2 ng/L, range UD to 10 ng/L
- #Algebraic Equations
- #pth
- pthp = K_PTH * (1.0/(1.0 + (ca/SetPoint_PTH_50))) - dpth * pth
- #dhcc
- dthp = K_DHCC*pth - ddhhc*dthh
- #ctn
- ctnp = K_CTN * (ca/SetPoint_CT) / (1.0 + (ca/SetPoint_CT) ) - dctn * ctn
- #calcium
- cap = K1*pth - K2*ctn + K3*dthh - ca*K4*((ctn/pth)/(1.0 + (ctn/pth))) - K5*ca
- n = len(x) # 4: implies we have 4 ODEs
- xp = np.zeros((n,1))
- xp[0] = cap
- xp[1] = pthp
- xp[2] = dthp
- xp[3] = ctnp
- return xp
- # The ''driver'' that will integrate the ODE(s):
- if __name__ == '__main__':
- # Start by specifying the integrator:
- # use ''vode'' with "backward differentiation forumla"
- r = integrate.ode(calcium).set_integrator('dopri5')
- #Set the time range
- t_start = 0.0
- t_final = 10.0
- delta_t = 0.001
- #Number of time steps: 1 extra for initial condition
- num_steps = np.floor((t_final - t_start)/delta_t) + 1
- #Set initial condition(s): for integrating variable and time!
- ca_zero = 95 #mg/L
- pth_zero = 35.8e-6 #mg/L
- dthh_zero = 2.163e-2 #mg/L
- ctn_zero = 2e-6 #mg/L
- r.set_initial_value([ca_zero,pth_zero,dthh_zero,ctn_zero], t_start)
- #Additional Python step: create vectors to store trajectories
- t = np.zeros((num_steps, 1))
- ca = np.zeros((num_steps, 1))
- pth = np.zeros((num_steps,1))
- dthh = np.zeros((num_steps,1))
- ctn = np.zeros((num_steps,1))
- t[0] = t_start
- ca[0] = ca_zero
- pth[0] = pth_zero
- dthh[0] = dthh_zero
- ctn[0] = ctn_zero
- #Integrate the ODE(s) across each delta_t timestep
- k = 1
- while r.successful() and k < num_steps:
- r.integrate(r.t + delta_t)
- #Store the results to plot later
- t[k] = r.t
- ca[k] = r.y[0]
- pth[k] = r.y[1]
- dthh[k] = r.y[2]
- ctn[k] = r.y[3]
- k += 1
- #All done! Plot the trajectories
- fig = figure()
- ax1 = subplot(411)
- ax1.plot(t, ca)
- ax1.legend(["ca"])
- ax1.grid('on')
- ax2 = subplot(412)
- ax2.plot(t,pth)
- ax2.grid('on')
- ax2.legend(["pth"])
- ax3 = subplot(413)
- ax3.plot(t,dthh)
- ax3.grid('on')
- ax3.legend(["dthh"])
- ax4 = subplot(414)
- ax4.plot(t,ctn)
- ax4.legend(["ctn"])
- ax4.grid('on')
- show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement