Dr. Maxime,
I have a message from Dr. Maxime below:
@Maxime:
I made that code a long time ago, which is in matlab:
maybe I 've modified some part of the equation, check if it correspond to the python code.
clear all
close all
dt=0.05;
tfinal = 20000;
t = 0 : dt : tfinal;
%param init
Gk = 0.064;
Gkca = 0.08;
Gl = 0.0055;
Kd = 0.01;
fc = 0.4;
alpha = 4 * 10^-5;
Kca = 0.01;
El = -20;
Ek = -83;
Gca2 = 0.02694061;
vca2 = -20.07451779;
Rca = 5.97139101;
Jbase = 0.02397327;
Cm = 1;
R = 8.314;
T = 310;
F = 96.487;
Ca0 = 3 * ones(length(t),1);
Istim = sin(23.1416(t-5000)/(tfinal))'1;%ones(length(t),1)0.1; %
%Istim =zeros(length(t),1);
%Istim(round(tfinal/dt0.1):round(tfinal/dt0.9)) = ones(round(tfinal/dt0.9)-round(tfinal/dt0.1)+1,1)0.25;
% resultat
Vm = zeros(length(t),1);
%Vm(1)=-50;
%Vm(1:round(1/dt))= ones(round(1/dt),1)0;
nk = zeros(length(t),1);
nk(1)=0.079257;
Ca = zeros(length(t),1);
Ca(1)=0.0001;
%force parameter hai murphy parametre
Ca05MLCK=.0004640758;%mM
nM=4.7135;
K2=0.139910^-3;%1/s
h=15.6;% nm
fp1 = 14.449610^-3;%1/s
gp1 = 3.612410^-3;%1/s
g1 = 0.13410^-3;%1/s
fpsigma =fp1; % sigma = h
gpsigma =gp1; % sigma = h
gsigma = g1 ; % sigma = h
K=5.0859;
K1 = zeros(length(t),1);
%M + Mp + AMp + AM = 1
M = ones(length(t),1);
Mp = zeros(length(t),1);
AMp = zeros(length(t),1);
AM = zeros(length(t),1);
Force = zeros(length(t),1);
for i = 1 : length(t)
Eca = ((R * T) / (2 * F)) * log(Ca0(i) / Ca(i));
hki = 1 / (1 + exp((4.2 - Vm(i)) / 21.1));
tnk = 23.75 * exp(-Vm(i) / 72.15);
Ica2 = Jbase + Gca2 * (Vm(i) - Eca) / (1 + exp(-(Vm(i) - vca2) / Rca));
Ik = Gk * nk(i) * (Vm(i) - Ek);
Ikca = Gkca * Ca(i)^2 / (Ca(i)^2 + Kd^2) * (Vm(i) - Ek);
Il = Gl * (Vm(i) - El);
Vm(i+1) = ((Istim(i) - Ica2 - Ik - Ikca - Il) / Cm)*dt + Vm(i);
nk(i+1) = ((hki - nk(i)) / tnk) * dt + nk(i);
Ca(i+1) = (fc * (-alpha * Ica2 - Kca * Ca(i)))*dt + Ca(i);
% force Hai Murphy
K1(i) = (Ca(i)*10^-1)^nM / (Ca05MLCK^nM + (Ca(i)*10^-1)^nM)*10^-3;
M(i+1) = (-K1(i) * M(i) + K2 * Mp(i) + gsigma * AM(i))*dt+M(i) ;
Mp(i+1) =( K1(i) * M(i) - (K2 + fpsigma)* Mp(i) + gpsigma * AMp(i))*dt+Mp(i);
AMp(i+1) = (K1(i) * AM(i) + fpsigma* Mp(i) - (gpsigma + K2) * AMp(i))*dt+AMp(i);
AM(i+1) = (-(K1(i)+gsigma) * AM(i) +K2 * AMp(i))*dt+AM(i);
Force(i+1) =(AMp(i+1) + AM(i+1)); %K *
if( mod(i,1000)==0)
i/length(t) % progress
end
end
figure(1)
subplot(411)
plot(t,Vm(1:end-1));title(‘Vm’);
subplot(412)
plot(t,nk(1:end-1));title(‘nk’);
subplot(413)
plot(t,Ca(1:end-1));title(‘Ca’);
subplot(414)
plot(t,Force(1:end-1));title(‘force’);
figure(2)
plot(Vm,Ca);title(‘Vm/Ca’)
figure(3)
subplot(311)
plot(t,K1);title(‘K1’)
subplot(312)
plot(t,M(1:end-1),‘r’);hold on;plot(t,Mp(1:end-1),‘b’);plot(t,AMp(1:end-1),‘g’);plot(t,AM(1:end-1),‘y’);
plot(t,M(1:end-1)+Mp(1:end-1)+AMp(1:end-1)+AM(1:end-1),‘c’)
title(‘M,Mp,AMp,AM’)
subplot(313)
plot(t,Force(1:end-1));title(‘F’)
You can remove the force model and jjsut take the electrical equations
Eca = ((R * T) / (2 * F)) * log(Ca0(i) / Ca(i));
hki = 1 / (1 + exp((4.2 - Vm(i)) / 21.1));
tnk = 23.75 * exp(-Vm(i) / 72.15);
Ica2 = Jbase + Gca2 * (Vm(i) - Eca) / (1 + exp(-(Vm(i) - vca2) / Rca));
Ik = Gk * nk(i) * (Vm(i) - Ek);
Ikca = Gkca * Ca(i)^2 / (Ca(i)^2 + Kd^2) * (Vm(i) - Ek);
Il = Gl * (Vm(i) - El);
Vm(i+1) = ((Istim(i) - Ica2 - Ik - Ikca - Il) / Cm)*dt + Vm(i);
nk(i+1) = ((hki - nk(i)) / tnk) * dt + nk(i);
Ca(i+1) = (fc * (-alpha * Ica2 - Kca * Ca(i)))*dt + Ca(i);
Dr Maxime Yochum
LTSI - UMR 1099 Inserm - Université de Rennes 1
SESAME (Epileptogenic Systems : Signal and Models)
Campus Baulieu - Bât. 22
35042 Rennes Cedex - France -
Email - maxime.yochum@univ-rennes1.fr
Thank you for time Dr. Francois,
Best regards
Saeed