Dear Dr. Francois,
Hope you are fine,
Is it possible to generate a simulation in Brainstorm by inserting the equations
of uterus electrical model?
It just a reduced model with 3 simple equations,
If yes how can I import these equations?
Best regards
Saeed
Have you looked at the information in these tutorials and forum posts?
http://neuroimage.usc.edu/brainstorm/Tutorials/TimeFrequency#Simulation
Dear Dr. Francois,
We have the equations of uterine activity in python like that:
def init(self,):
self.Gk = 0.064
self.Gkca = 0.08
self.Gl = 0.0055
self.Kd = 0.01
self.fc = 0.4
self.alpha = 4 * 10 ** -5
self.Kca = 0.01
self.El = -20
self.Ek = -83
self.Gca2 = 0.02694061
self.vca2 = -20.07451779
self.Rca = 5.97139101
self.Jbase = 0.02397327
self.R = 8.314
self.T = 295
self.F = 96.487
self.Ca0 = 3
self.Istim = 0.
self.cm = 1
def derivT(self, dt, Vm, nk, Ca):
"""Computes temporal derivative for *Red3* model."""
#Nerst
Eca = ((self.R * self.T) / (2 * self.F)) * numpy.log(self.Ca0 / Ca)
#H inf x
hki = 1 / (1 + numpy.exp((4.2 - Vm) / 21.1))
#Tau x
tnk = 23.75 * numpy.exp(-Vm / 72.15)
#Courants
Ica2 = self.Jbase + self.Gca2 * (Vm - Eca) / \
(1 + numpy.exp(-(Vm - self.vca2) / self.Rca))
Ik = self.Gk * nk * (Vm - self.Ek)
Ikca = self.Gkca * Ca ** 2 / (Ca ** 2 + self.Kd ** 2) * (Vm - self.Ek)
Il = self.Gl * (Vm - self.El)
#Derivees
Vm = ((self.Istim - Ica2 - Ik - Ikca - Il) / self.cm) * dt + Vm
nk = ((hki - nk) / tnk) * dt + nk
Ca = (self.fc * (-self.alpha * Ica2 - self.Kca * Ca)) * dt + Ca
return Vm, nk, Ca
Thank you for suggestions,
Best regards
Saeed
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
I’m sorry, I cannot help you with this.
This is out of my range of competence.
OK Dr. Francois no problem but result of the code should be like the attached code
this Vm for one node.