calciumregPython

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 = 1.6e6 #PTH
    K2 = -18e6 #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()