# This Python file describes a collection of electrically coupled
# pseudo-plateau bursting cells as described in the paper
# "Network Properties of Electrically Coupled Bursting Pituitary
# Cells", Mehran Fazli and Richard Bertram, Frontiers in Endocrinology,
# 13:936160, 2022.
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
import seaborn as sns
from random import random
def ode_solver(IC,Nt,ggj):
# parameter values
vn=-5
kc=0.12
ff=0.005
vca=60
vk=-75
vl = -50.0
gk=2.5
cm=5
gbk=1
gca=2.1
gsk=2
vm=-20
vb=-5
sn=10
sm=12
sbk=2
taun=30
taubk=5
ks=0.4
alpha=0.0015
gl=0.2
# ggj=0.002
dt=0.5
#//////////////////////////////////
# LOOP OVER INITIAL CONDITIONS
vtim=np.zeros((num_cell,Nt))
ntim=np.zeros((num_cell,Nt))
ctim=np.zeros((num_cell,Nt))
btim=np.zeros((num_cell,Nt))
vv=np.zeros(num_cell)
nn=np.zeros(num_cell)
cc=np.zeros(num_cell)
bb=np.zeros(num_cell)
T_n=np.zeros(num_cell)
T_ij=np.zeros((num_cell,num_cell))
C_ij=np.zeros((num_cell,num_cell))
#/////////////// initialize the variable matrices
for i in range(0,num_cell):
vv[i]= IC[0,i]
nn[i]= IC[1,i]
cc[i]= IC[2,i]
bb[i]= IC[3,i]
for t in range(0,Nt): #/////////////////////////////////////////////////SOLVER TIME-STEP ITERATION
ninf=0
bkinf=0
minf=0
cinf=0
ica=0
isk=0
ibk=0
ikdr=0
il=0
k1v=0
k2v=0
k3v=0
k4v=0
k1n=0
k2n=0
k3n=0
k4n=0
k1c=0
k2c=0
k3c=0
k4c=0
k1b=0
k2b=0
k3b=0
k4b=0
vv_old=np.zeros(num_cell)
for i in range(0,num_cell):
vv_old[i]=vv[i]
for i in range(0,num_cell): #/////////////////////// LOOP IN CELLS
#//////////////////////////////////////////////////////////////////////////// FIRST STEP RK
ninf=1/(1+np.exp((vn-vv[i])/sn))
bkinf=1/(1+np.exp((vb-vv[i])/sbk))
minf=1/(1+np.exp((vm-vv[i])/sm))
cinf=((cc[i])**2)/(((cc[i])**2)+ks*ks)
ica=gca*minf*(vv[i]-vca)
isk=gsk*cinf*(vv[i]-vk)
ibk=gbk*bb[i]*(vv[i]-vk)
ikdr=gk*nn[i]*(vv[i]-vk)
il = gl*(vv[i]-vl)
igj=0;
for h in range(0,num_cell):
igj=igj+SN_sim[i][h]*ggj*(vv[i]-vv_old[h])
k1v = dt*(-(ica+isk+ibk+ikdr+igj+il)/cm)
k1n = dt*((ninf-nn[i])/taun)
k1c = dt*(-ff*(alpha*ica+kc*cc[i]))
k1b = dt*(-(bb[i]-bkinf)/taubk)
#//////////////////////////////////////////////////////////////////////////// SECOND STEP RK
ninf=1/(1+np.exp((vn-(vv[i] + 0.5*k1v))/sn))
bkinf=1/(1+np.exp((vb-(vv[i] + 0.5*k1v))/sbk))
minf=1/(1+np.exp((vm-(vv[i] + 0.5*k1v))/sm))
cinf=((cc[i] + 0.5*k1c)**2)/(((cc[i] + 0.5*k1c)**2)+ks*ks)
ica=gca*minf*((vv[i] + 0.5*k1v)-vca)
isk=gsk*cinf*((vv[i] + 0.5*k1v)-vk)
ibk=gbk*(bb[i] + 0.5*k1b)*((vv[i] + 0.5*k1v)-vk)
ikdr=gk*(nn[i] + 0.5*k1n)*((vv[i] + 0.5*k1v)-vk)
il = gl*((vv[i] + 0.5*k1v)-vl)
igj=0
for h in range(0,num_cell):
igj=igj+SN_sim[i][h]*ggj*((vv[i] + 0.5*k1v)-vv_old[h])
k2v = dt*(-(ica+isk+ibk+ikdr+igj+il)/cm)
k2n = dt*((ninf-(nn[i] + 0.5*k1n))/taun)
k2c = dt*(-ff*(alpha*ica+kc*(cc[i] + 0.5*k1c)))
k2b = dt*(-((bb[i] + 0.5*k1b)-bkinf)/taubk)
#//////////////////////////////////////////////////////////////////////////////// THIRD STEP RK
ninf=1/(1+np.exp((vn-(vv[i] + 0.5*k2v))/sn))
bkinf=1/(1+np.exp((vb-(vv[i] + 0.5*k2v))/sbk))
minf=1/(1+np.exp((vm-(vv[i] + 0.5*k2v))/sm))
cinf=((cc[i] + 0.5*k2c)**2)/(((cc[i] + 0.5*k2c)**2)+ks*ks)
ica=gca*minf*((vv[i] + 0.5*k2v)-vca)
isk=gsk*cinf*((vv[i] + 0.5*k2v)-vk)
ibk=gbk*(bb[i] + 0.5*k2b)*((vv[i] + 0.5*k2v)-vk)
ikdr=gk*(nn[i] + 0.5*k2n)*((vv[i] + 0.5*k2v)-vk)
il = gl*((vv[i] + 0.5*k2v)-vl);
igj=0;
for h in range(0,num_cell):
igj=igj+SN_sim[i][h]*ggj*((vv[i] + 0.5*k2v)-vv_old[h])
k3v = dt*(-(ica+isk+ibk+ikdr+igj+il)/cm)
k3n = dt*((ninf-(nn[i] + 0.5*k2n))/taun)
k3c = dt*(-ff*(alpha*ica+kc*(cc[i] + 0.5*k2c)))
k3b = dt*(-((bb[i] + 0.5*k2b)-bkinf)/taubk)
#////////////////////////////////////////////////////////////////////////////////// FOURTH STEP RK
ninf=1/(1+np.exp((vn-(vv[i] + 0.5*k3v))/sn))
bkinf=1/(1+np.exp((vb-(vv[i] + 0.5*k3v))/sbk))
minf=1/(1+np.exp((vm-(vv[i] + 0.5*k3v))/sm))
cinf=((cc[i] + 0.5*k3c)**2)/(((cc[i] + 0.5*k3c)**2)+ks*ks)
ica=gca*minf*((vv[i] + 0.5*k3v)-vca)
isk=gsk*cinf*((vv[i] + 0.5*k3v)-vk)
ibk=gbk*(bb[i] + 0.5*k3b)*((vv[i] + 0.5*k3v)-vk)
ikdr=gk*(nn[i] + 0.5*k3n)*((vv[i] + 0.5*k3v)-vk)
il = gl*((vv[i] + 0.5*k3v)-vl)
igj=0;
for h in range(0,num_cell):
igj=igj+SN_sim[i][h]*ggj*((vv[i] + 0.5*k3v)-vv_old[h])
k4v = dt*(-(ica+isk+ibk+ikdr+igj+il)/cm)
k4n = dt*((ninf-(nn[i] + 0.5*k3n))/taun)
k4c = dt*(-ff*(alpha*ica+kc*(cc[i] + 0.5*k3c)))
k4b = dt*(-((bb[i] + 0.5*k3b)-bkinf)/taubk)
#////////////////////////////////////////////////////////////////////////////////// FINAL STEP RK
vv[i] = vv[i] + (1.0/6.0)*(k1v + 2*k2v + 2*k3v + k4v)
nn[i] = nn[i] + (1.0/6.0)*(k1n + 2*k2n + 2*k3n + k4n)
cc[i] = cc[i] + (1.0/6.0)*(k1c + 2*k2c + 2*k3c + k4c)
bb[i] = bb[i] + (1.0/6.0)*(k1b + 2*k2b + 2*k3b + k4b)
vtim[i,t]=vv[i]
ntim[i,t]=nn[i]
ctim[i,t]=cc[i]
btim[i,t]=bb[i]
if t>=Nt-20000:
for i in range(num_cell-1):
for h in range(i+1,num_cell):
if vv[i]>=-40 and vv[h]>=-40:
T_ij[i,h]=T_ij[i,h]+1
if i==0 and h==1:
if vv[i]>=-40:
T_n[i]=T_n[i]+1
if vv[h]>=-40:
T_n[h]=T_n[h]+1;
if i==0 and h>1:
if vv[h]>=-40:
T_n[h]=T_n[h]+1
for i in range(num_cell-1):
for h in range(i+1,num_cell):
C_ij[i,h]=T_ij[i,h]/((T_n[i]*T_n[h])**0.5)
C_ij[h,i]=T_ij[i,h]/((T_n[i]*T_n[h])**0.5)
return C_ij, vtim, ntim, ctim, btim
# a two node network
SN_sim=[[0, 1],
[1, 0]];
G_sim_sc = nx.Graph()
num_cell=len(SN_sim)
for i in range(num_cell):
for j in range(num_cell):
if SN_sim[i,j] == 1:
G_sim_sc.add_edge(i,j)
ggj=0.002 # gap junction conductance
dt=0.5
sec=100
num_set=100 # number of initial conditions
Ntim=int((sec*1000*(1/dt)))
sim=np.zeros((num_set,num_cell,num_cell))
ICall=np.zeros((num_set,4,num_cell))
vall=np.zeros((num_set,num_cell,Ntim))
nall=np.zeros((num_set,num_cell,Ntim))
call=np.zeros((num_set,num_cell,Ntim))
ball=np.zeros((num_set,num_cell,Ntim))
for setn in range(1,num_set+1):
for i in range(0,num_cell):
ICall[setn-1,0,i]= -90*(random())+20
ICall[setn-1,1,i]= random()/2
ICall[setn-1,2,i]= random()/2
ICall[setn-1,3,i]= random()
AA=ICall[setn-1,:,:]
sim[setn-1,:,:], vall[setn-1,:,:], nall[setn-1,:,:], call[setn-1,:,:], ball[setn-1,:,:]=ode_solver(AA,Ntim,ggj)