- Author:
- aram148 <42922407+aram148@users.noreply.github.com>
- Date:
- 2022-07-22 15:47:05+12:00
- Desc:
- Added documentation for VSM model
- Permanent Source URI:
- https://models.cellml.org/workspace/6b0/rawfile/a8a92308e217ac5626809237dd90a31240b22834/USMC/DM_USM.py
import numpy as np
import scipy as sp
import math
import matplotlib.pyplot as plt
from numpy import loadtxt
import opencor as oc
simulation = oc.open_simulation(r"C:\Users\aram148\Desktop\Platform1\EP3\USMC\Vm_stim_experiment.cellml")
# data = simulation.data()
# data.set_starting_point(0)
# data.set_ending_point(120)
# dt = (120-0)/10000
# data.set_point_interval(dt)
# simulation.run()
# results=simulation.results()
# k11 = results.algebraic()['outputs/K_1'].values()
# k1 = np.random.random((len(k11)-1,1))
# k1 = loadtxt('k1.csv', delimiter=',')
def DM_funcs_USM(t,k1,R):
# def DM_funcs_USM(t,R):
# k1 = 5.166796617067728e-05
# k2 = 1.2387
# g1 = 0.0756
# gp1 = 0.0709
# fp1 = 0.2838
# k1 = 0.9;
k2 = 0.1399
fp1 = 14.4496
gp1 = 3.6124
g1 = 0.1340
r = R
gam = 0
# print(r[0])
# Mean for cdf
p1 = r[1]/r[0]
p2 = r[4]/r[3]
# Std dev for cdf
q1 = np.sqrt((r[2]/r[0])-(r[1]/r[0])**2)
q2 = np.sqrt((r[5]/r[3])-(r[4]/r[3])**2)
# r, phi and I for M1_lambda
r0 = -p1/q1
r1 = (1-p1)/q1
phi0 = 0.5*(1+math.erf((r0-p1)/(q1*np.sqrt(2))))
phi1 = 0.5*(1+math.erf((r1-p1)/(q1*np.sqrt(2))))
I0 = -(np.exp(-((-p1/q1)**2)/2))/(np.sqrt(2*np.pi))
I1 = -(np.exp(-((((1-p1)/q1))**2)/2))/(np.sqrt(2*np.pi));
# r, phi and I for M2_lambda
r20 = -p2/q2
r21 = (1-p2)/q2
phi20 = 0.5*(1+math.erf((r20-p2)/(q2*np.sqrt(2))))
phi21 = 0.5*(1+math.erf((r21-p2)/(q2*np.sqrt(2))))
I20 = -(np.exp(-((-p2/q2)**2)/2))/(np.sqrt(2*np.pi))
I21 = -(np.exp(-((((1-p2)/q2))**2)/2))/(np.sqrt(2*np.pi))
# Functions for the RHS of M1 PDE
J0 = phi0
J10 = ((p1*phi0) + (q1*I0))
J11 = (p1*phi1) + (q1*I1)
J20=((p1**(2))*phi0)+((2*p1*q1)*I0)+((q1**(2))*(phi0+(r0*I0)))
J21=((p1**(2))*phi1)+((2*p1*q1)*I1)+((q1**(2))*(phi1+(r1*I1)))
J30=(p1**(3)*phi0)+((3*p1**(2)*q1)*I0)+((3*p1*q1**(2))*((phi0)+(r0*I0)))+((q1**(3)*(2+r0**(2))*I0))
J31=(p1**(3)*phi1)+((3*p1**(2)*q1)*I1)+((3*p1*q1**(2))*(phi1+(r1*I1)))+(q1**(3)*(2+(r1**(2)))*I1)
#Functions defined for the RHS of the second PDE M2_lambda
K0 = phi20
K01 = phi21
K10 = (p2*phi20) + (q2*I20)
K11 = (p2*phi21) +( q2*I21)
K20=((p2**(2))*phi20) + ((2*p2*q2)*I20) + ((q2**(2))*(phi20 + (r20*I20)))
K21=((p2**(2))*phi21) + ((2*p2*q2)*I21) + ((q2**(2))*(phi21 + (r21*I21)))
K30=(p2**(3)*phi20) + ((3*p2**(2)*q2)*I20) + ((3*p2*q2**(2))*((phi20) + (r20*I20))) + ((q2**(3)*(2+r20**(2))*I20))
K31=(p2**(3)*phi21) + ((3*p2**(2)*q2)*I21) + ((3*p2*q2**(2))*(phi21 + (r21*I21))) + (q2**(3)*(2+(r21**(2)))*I21)
# Components for the matrix F that will represent each moment,
# M1_lambda and M2_lambda
A0 = ((fp1*(1-r[6]))/2)-(fp1*(J11-J10)*r[0])
A1 = ((fp1*(1-r[6]))/3)-(fp1*(J21-J20)*r[0])
A2 = ((fp1*(1-r[6]))/4)-(fp1*(J31-J30)*r[0])
B0=(3*(fp1+gp1)*J0)+gp1*(J11-J10)+(4*gp1)*(p1-J11)
B1=(3*(fp1+gp1)*J10)+gp1*(J21-J20)+(4*gp1)*((p1**(2)+q1**(2))-J21)
B2=(3*(fp1+gp1)*J20)+gp1*(J31-J30)+(4*gp1)*((p1**(3)+3*p1*q1**(2))-J31)
C0=k1
C1=k1*p2
C2=k1*(p2**(2)+q2**(2))
D0=k2
D1=k2*p1
D2=k2*(p1**(2)+q1**(2))
E0=(20*g1*K0)+g1*(K11-K10)+(g1)*(1-K01)
E1=(20*g1*K10)+g1*(K21-K20)+(g1)*((p2)-K11)
E2=(20*g1*K20)+g1*(K31-K30)+(g1)*(p2**(2)+q2**(2)-K21)
V = gam*(A1-E0-B0)/(1+gam*(r[0]+r[3]))
F=np.array([A0-B0*r[0]+C0*r[3]-k2*r[0], A1-B1*r[0]+C1*r[3]-k2*r[1]-V*r[0], A2-B2*r[0]+C2*r[3]-k2*r[2]-2*V*r[1], D0*r[0]-E0*r[3]-k1*r[3], D1*r[0]-E1*r[3]-k1*r[4]-V*r[3], D2*r[0]-E2*r[3]-k1*r[5]-2*V*r[4], -k1*r[6]+(1-r[6])*k2])
# simulation.clear_results()
return F
def USM_DM(trange,N,M,xrange):
T = np.linspace(trange[0], trange[1], N)
X0 = np.linspace(xrange[0], xrange[1], M)
dt=(trange[1]-trange[0])/(N)
M10=0.005
M11=0.001
M12=0.01
M20=0.005
M21=0.001
M22=0.01
C=1
# M10=0.5140
# M11=0.3110
# M12= 0.2146
# M20=0.1441
# M21=0.0870
# M22=0.0597
# C=0.2222
# M10 = 0.341238863408273
# M11 = 0.243924139636747
# M12 = 0.174361692884738
# M20 = 0.0910630566362484
# M21 = 0.0620757257099705
# M22 = 0.0478612210259969
# C = 0.352924011750745
# R = np.zeros(7)
# Calculating the Ca2+ for the same trange
simulation.reset(True)
data = simulation.data()
data.set_starting_point(trange[0])
data.set_ending_point(trange[1])
dt = (trange[1]-trange[0])/N
data.set_point_interval(dt)
simulation.run()
results=simulation.results()
k1 = results.algebraic()['outputs/K_1'].values()
Cai = results.states()['outputs/Cai'].values()
Vm = results.algebraic()['mPulse_protocol_s/V'].values()
R = np.array([M10, M11, M12, M20, M21, M22, C])
Rnew = np.zeros([len(T),len(R)])
Rnew[0,:] = R
# i = 0
for i in range(1,len(T)):
t = T[i-1]
trange2 = T[i]
F1=DM_funcs_USM(t,k1[i-1],Rnew[i-1,:])
F2 =DM_funcs_USM(t+(dt/2),k1[i-1],Rnew[i-1,:]+(dt/2)*F1)
F3 =DM_funcs_USM(t+(dt/2),k1[i-1],Rnew[i-1,:]+(dt/2)*F2)
F4 =DM_funcs_USM(t+(dt),k1[i-1],Rnew[i-1,:]+(dt)*F3)
# F1 =DM_funcs_USM(t,Rnew[i-1,:])
# F2 =DM_funcs_USM(t+(dt/2),Rnew[i-1,:]+(dt/2)*F1)
# F3 =DM_funcs_USM(t+(dt/2),Rnew[i-1,:]+(dt/2)*F2)
# F4 =DM_funcs_USM(t+(dt),Rnew[i-1,:]+(dt)*F3)
Rnew[i,:] = Rnew[i-1,:]+ (dt/6)*(F1 + (2*F2) + (2*F3) + F4)
# print(i)
# print(Rnew[i,0])
# Rnew[i,:] = Rnew[i-1,:] + dt*F1
return Rnew, Cai, Vm, k1
# F = DM_funcs_USM(t, k1, R)
# RR = lambda t, s:F
# sol =solve_ivp(RR,T,R,method='RK4')
# return sol
trange = [0, 20]
xrange = [-20, 20]
N = 50000
M = 1000
kappa = 5.0859
gam =np.array([0, 44.87, 141.03])
T = np.linspace(trange[0], trange[1], N)
Rnew, Cai, Vm, k1 = USM_DM(trange, N, M, xrange)
force = (Rnew[:,1]+Rnew[:,4])
stiffness = Rnew[:,0]+Rnew[:,3]
nmp = 1-Rnew[:,6]-Rnew[:,0]
phosphorylation = nmp + Rnew[:,0]
figure, axis = plt.subplots(3, 2)
axis[0,0].plot(T,force,'red')
axis[0, 0].set_title("Force")
axis[0,1].plot(T,stiffness,'blue')
axis[0, 1].set_title("stiffness")
axis[1,0].plot(T,phosphorylation,'green')
axis[1, 0].set_title("phosphorylation")
axis[1,1].plot(T,k1[1:],'yellow')
axis[1, 1].set_title("k1")
axis[2,0].plot(T,Cai[1:])
axis[2, 0].set_title("Ca_i")
axis[2,1].plot(T,Vm[1:])
axis[2,1].set_title("Vm")
axis[1,1].set_title("k1")
# plt.plot(T,force)
plt.show()
