The model equations are expressed as follows (Gerstner and
Kistler, 2002):
C
m
dV
i
=dt ¼g
Na
m
3
i
h
i
ðV
i
V
Na
Þg
K
n
4
i
ðV
i
V
K
Þg
L
ðV
i
V
L
Þ
þI
ext
i
þ I
syn
i
þ I
noise
i
þ γ
i
I
astro
i
;
dn
i
=dt ¼ α
n
i
ðV
i
Þð1n
i
Þβ
n
i
ðV
i
Þn
i
;
dm
i
=dt ¼ α
m
i
ðV
i
Þð1m
i
Þβ
m
i
ðV
i
Þm
i
;
dh
i
=dt ¼ α
h
i
ðV
i
Þð1h
i
Þβ
h
i
ðV
i
Þh
i
; ð1Þ
where i ¼1; 2; 3 denote the neuron index, V
i
represent the mem-
brane potentials for each neuron i in the motif, g
K
, g
Na
, and g
L
stand for the maximal potassium, sodium, and leakage conduc-
tances per unit area, and V
K
, V
Na
, and V
L
denote their correspond-
ing reversal potentials. Without loss of generality, the maximum
conductance of the potassium, solidum, and leakage channels are
set as g
K
¼36 ms/cm
2
, g
Na
¼120 ms/cm
2
and g
L
¼0:3 ms/cm
2
,
respectively, and their corresponding reversal potentials are set as
V
K
¼12 mV, V
Na
¼115 mV and V
L
¼10:6 mV, respectively. The
membrane capacitance is set as C
m
¼1 μ F/cm
2
. The parameters m
i
,
h
i
, n
i
represent the saturation values of the gating variables with
those terms in (1) governed by
α
n
i
¼0:01ðV
i
þ 50Þ=½1expð0:1ðV
i
þ 50ÞÞ;
β
n
i
¼0:125 expð0:125ðV
i
þ 60ÞÞ;
α
m
i
¼0:1ðV
i
þ 35Þ=½1expð0:1ðV
i
þ 35ÞÞ;
β
m
i
¼4:0 expð0:0556ðV
i
þ 60ÞÞ;
α
h
i
¼0:07 expð0:05ðV
i
þ 60ÞÞ;
β
h
i
¼1=½1 þ expð0:1ðV
i
þ 30ÞÞ:
In Eq. (1), the term I
ext
i
denotes the externally applied current.
The term I
syn
i
denotes the synaptic coupling current among the
neurons, which is modeled by the linear sum of the synaptic
current onto neuron i coming from all the pre-synaptic neurons,
i.e. I
syn
i
¼∑
j
g
ij
r
j
ðEV
i
Þ, where g
ij
describes the coupling strength of
the synapse from neuron j to neuron i. For simplicity, we assume
the coupling strength is identical for all the connections, denoted
as g
ij
¼g for short. The reversal potential E¼0 mV if neuron j is
excitatory and E ¼80 mV if neuron j is inhibitory. When a
neuron fires an action potential, the synaptic conductance r
j
of its postsynaptic target neuron is increased according to
r
j
¼r
j
þ 0:5ð1r
j
Þ, and in other time dr
j
=dt ¼r
j
=τ
2
, where the
refractory time τ
2
is set as 15 ms. The term I
noise
i
denotes the noisy
synaptic current representing the external or internal fluctuations,
which is modeled by the Ornstein–Uhlenbeck (OU) process
(Destexhe and Rudolph, 2012)
τ
d
dI
noise
i
=dt ¼I
noise
i
þ
ffiffiffiffiffiffiffi
2D
p
ηðtÞ; ð2Þ
where ηðtÞ is a Gaussian white noise with zero mean and unit
variance. For simplicity, here we assume the noise intensity for
each neuron is the same, denoted as D. The constant τ
d
denotes the
relaxation time of the synaptic noise, which is fixed at τ
d
¼2 ms.
The synaptic inward current I
astro
i
in (1) denotes the depolar-
ization current feedbacked from astrocyte to neuron i with a
coupling strength γ
i
. The value of I
astro
i
is related to the calcium
concentration in a two-dimensional field. More details are given in
Section 2.2.
2.2. Astrocyte field model
In literature, several models have been suggested to describe
the dynamics of the intracellular Ca
2þ
concentration (Li and Rinzel,
1994; Atri et al., 1993; Amiria et al., 2012; Höfer et al., 2002; Sneyd
et al., 1995). In the study of dressed neuron (a neuron coupled
by an astrocyte), the Li–Rinzel model is typically used, which
resembles the Hodgkin–Huxley mode for electrically excitable
membranes by replacing the transmembrane potential by Ca
2þ
concentration. In this model, the calcium only changes in the time
domain, and the signaling way from astrocyte to neuron is directly
modeled by a depolarization current determined by the time-
dependent Ca
2þ
.InSotero and Cancino (2010), a dynamical mean
field was proposed to model the behavior of an ensemble of
interacting neurons and astrocytes, which gives a simple and
efficient description of neural-glial mass. Yet, there still exist
rooms for modeling improvement. Firstly, the astrocytes are
spatially extended objects. The Ca
2þ
can be diffused spatially from
astrocyte to astrocyte through gap junctions, and thus the Ca
2þ
concentration in different spatial domains may be different
(Allegrini et al., 2009; Liu and Li, 2013). The original Li–Rinzel
model is unable to represent the Ca
2þ
concentration in spatial
domain and thus cannot give a more accurate evaluation of the
depolarization current I
astro
i
, as it is position-dependent. Secondly,
for a network motif coupled by several astrocytes, it is difficult to
model the precise relations between astrocytes and neurons if
only the time domain is considered, since each astrocyte may
interact with several neurons in its neighborhood and meanwhile
each neuron may be affected by more than one astrocyte
(Newman, 2003). Considering these two aspects, it is more
reasonable to model the astrocytes as a two-dimensional field,
and then the depolarization current from astrocyte to neuron can
be directly determined by Ca
2þ
in the spatial domain (Allegrini
et al., 2009; Liu and Li, 2013).
Following this motivation, in this paper, the Li–Rinzel model is
extended to a two-dimensional field, and the diffusion effects of
gliotransmitters, IP
3
and Ca
2þ
from astrocyte to astrocyte are
introduced. Mathematically, the model equations can be expressed
as follows (Ullah et al., 2006):
∂½IP
3
=∂t ¼D
p
ð∂
2
½IP
3
=∂
2
x þ ∂
2
½IP
3
=∂
2
yÞþ1=τ
p
ð½IP
3
n
½IP
3
Þ
þ½v
p
ð½Ca
2þ
þð1αÞk
p
Þ=ðk
p
þ½Ca
2þ
Þ þ r
j
ΘðV
j
25mVÞ
ð3aÞ
∂½Ca
2þ
=∂t ¼ D
c
ð∂
2
½Ca
2þ
=∂
2
x þ ∂
2
½Ca
2þ
=∂
2
yÞJ
flux
J
leak
J
pump
ð3bÞ
τ
q
dq=dt ¼ k
2
2
=ðk
2
2
þ½Ca
2þ
2
Þq: ð3cÞ
The first terms in right hand side (RHS) of (3a) and (3b) reflect
the diffusion of IP
3
and Ca
2þ
, moving from astrocyte to astrocyte
via gap junction with diffusion coefficients D
p
¼300ðμmÞ
2
s
1
and
D
c
¼20ðμmÞ
2
s
1
, respectively (Allegrini et al., 2009; Atri et al.,
1993). In the dynamical equation (3a) for IP
3
, the second term
describes the degradation of IP
3
with a degradation rate of 1=τ
p
and an equilibrium point ð½IP
3
Þ
n
, which accounts for a refractory
time in calcium responses of astrocytes; the third term describes
the production of IP
3
caused by Ca
2þ
, where the constants are
v
p
¼0:13 μMs
1
, k
p
¼1:1 μM, and α ¼0:8; the last term reflects the
glutamate-induced production of the cytosolic IP
3
with the pro-
duction rate r
j
via the Heaviside function Θ. To be specific, it is
activated when the membrane potential of the neuron j,
V
j
≥25 mV. It is also worth pointing out that as the inhibitory
neurons cannot release much excitatory neurotransmitter, gluta-
mate, which is necessary for establishing neuron–astrocyte inter-
action, the signaling way from neuron to astrocyte is not
considered for the inhibitory neurons, i.e. r
j
¼0 if neuron j is
inhibitory.
In the dynamical equation (3b) for Ca
2þ
concentration, J
flux
describes the rate of calcium released from ER into the cytoplasm
through the IP
3
receptors, which designates a positive feedback
controlled by the slow inactivation gate q, J
pump
denotes the rate of
calcium flux pumping from the cytoplasm into ER and into the
extracellular space; J
leak
denotes the leakage flux from the ER to
Y. Liu, C. Li / Journal of Theoretical Biology 335 (2013) 265–275 267