# Simulation with new equations

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/dt
0.1):round(tfinal/dt
0.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.1399
10^-3;%1/s
h=15.6;% nm
fp1 = 14.4496
10^-3;%1/s
gp1 = 3.612410^-3;%1/s
g1 = 0.134
10^-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