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()