A PHYSICS - INFORMED NEURAL NETWORK FOR COUPLEDCALCIUM DYNAMICS IN A CABLE NEURON

Zachary M. Miksis* and Gillian Queisser†
Department of Mathematics
Temple University
Philadelphia, PA 19122

A BSTRACTTranscranial magnetic stimulation (TMS) is a noninvasive treatment for a variety of neurologicaland neuropsychiatric disorders by triggering a calcium response through magnetic stimulation. To
understand the full effects of this treatment, researchers will often use numerical simulations to modeland study the calcium response. These simulations are limited to short-time simulations of single
neurons due to computational complexity, restricting their use in clinical settings. In this paper, we
explore an application of physics-informed neural networks (PINNs) to accurately produce long-time
simulations of neuronal responses, opening the possibility of utilizing these methods in clinicalapplications to directly benefit patients.

1

Introduction

Transcranial magnetic stimulation (TMS) is an invaluable
tool for treating a variety of neurological and neuropsychiatric disorders in a noninvasive way by using a timevarying magnetic field passed through the brain to stimulate neurons (2; 7). This stimulation triggers a response
from intracellular calcium, which is vital to regulating thetransfer of information from synaptic sites to the cell nucleus (4). Ultimately, this allows for treatment of various
forms of degenerative neurological disorders, and TMS
is therefore used extensively in both research and clinicalsettings (8).
While TMS is very important, its full effects are still notentirely understood. To investigate these effects, computational tools and studies are exceptionally important tocomplement experimental studies (5). To this end, different numerical computation tools have been developed(3; 13; 6). While these tools are very good at capturing
fine grain details and the quick responses to fast TMS frequencies, they are limited to short simulated time spans of
single neurons due to computational complexity. In orderto address these limitations, we look to neural networks(from here on, "neuron" will be used to refer to a biologicalneuron, and "neural network" will be used to refer to themachine learning paradigm).
∗†

Neural networks have been used in various applications,many focusing on image recognition and reconstruction(e.g., see (1)). More recently, knowledge of physical models have been introduced into neural networks to form
physics-informed neural networks (PINNs) (11; 9; 10). In
this work, we utilize these PINNs to incorporate partial differential equations (PDEs) that model calcium dynamics inneurons into the larger neural network. With a network thatsatisfies physical laws, we can build simulations that provide long-time results, overcoming one of the limitationsof traditional numerical methods. This leads to anotherlimitation of PINNs, in that they have been shown to havedifficulty accurately to simulate diffusion problems (12).
This can be overcome with strategies of both relaxation of
the loss function and periodic activation functions (15; 14),which we utilize in our application.
This paper is structured as follows: In Section 2, we describe the physical model of calcium dynamics that is used(with further details in Appendix A). In Section 3, we
describe the details on PINNs as we used them in our application. In Section 4, we present results of a simulated
cable neuron. In Section 5, we draw conclusions from ourresults and describe the broader impacts of this work.

miksis@temple.edu
queisser@temple.edu

International Conference on Scientific Computing and Machine Learning 2024 (SCML2024)

Coupled Calcium PINN

2

Calcium Dynamics

with x ∈ Ω, t ∈ [0, T ], where u(x, t) is the latent solution,
N [·] is a nonlinear differential operator, and Ω is the computational domain. Define f (x, t) as the left-hand-side ofEq. (8),f := ut + N [u].
(9)
By approximating u(x, t) using a deep neural network,
we build a physics-informed neural network f (x, t) by
combining this network with Eq. (9). The shared terms ofthe networks are learned by minimizing the mean squared
error loss
M SE = (1 − α)M SEu + α(M SEf + M SEb ), (10)
where
Nu
1 X
M SEu =|u(tiu , xiu ) − ui |2(11)
Nu i=1

Within a neuron, ion dynamics are modeled by reactiondiffusion equations of the form∂u
= D∆u + R(u),∂t(1)

where u is an ion concentration, D is the diffusion coefficient, and R(u) is a reaction term. From Fick’s first law,
we have∂u
J = −D ,(2)∂xwhere J is the diffusion flux. This leads to the Neumannboundary condition
∂u
J
=− .
∂x
D

(3) represents the error of the initial condition,
Nf
1 X
M SEf =|f (tif , xif )|2(12)
In the cytosol, the non-organelle interior of a neuron,
Nf i=1
calcium is transported via various mechanisms and isbuffered by the molecule calbindin. We model this us- represents the error of the operator f inside the domain,anding a dimension-reduced system (4), given byNb
1 X|ux (tib , xib ) + J/D|2
(13)
M
SE
=
b
∂c
N
b
= ∇ · (Dc ∇cc ) + f (b, cc ) + JP M ,
(4)
i=1
∂t
represents the error of the Neumann boundary condi∂bu= ∇ · (Db ∇b) + f (b, cc ),(5) tions. Here, {tiu , xiu , ui }N
i=1 represents the initial condi∂tii Nftion training data, {tf , xf }i=1 represents the collocation
where the reaction term
points within the domain (excluding the boundaries), and
b
− tot
+
{tib , xib }N
f (b, cc ) = kb (b − b) − kb bcc
i=1 represents the boundary collocation points.
In our minimization, we apply a relaxation parametermodels the reaction equationα ∈ (0, 1]. This is selected a priori by trial and error,
and may require tuning to fit a particular model (15).
−
κb
2+
2+
−−
⇀
(6) A defining feature of the calcium model is that it is a
Ca + CalB ↽−− CalBCa ,κ+bpure diffusion problem, which have been difficult to accurately model with neural networks. PINNs were initally
2+
and JP M represents the net Ca
ion flux across the introduced with tanh activation function, and were demonplasma membrane (separating intra- and extracellular strated to perform well for advection-dominated problems
space). This term can be broken down into the differ- (11). The SIREN network introduced the use of periodic
ent mechanisms that transport calcium across the plasma activation functions (14), which have been demonstrated
membrane,to correct for problems in modeling diffusion-dominated
JP M = −JP − JN + JSY N + JV DCC ,(7) problems (12). Therefore, we utilize a sin activation function in our network to provide consistent long-time modelwhere JP is the flux from plasma membrane Ca2+ -ATPase ing.
pumps (PMCA), JN is the flux from Na+ /Ca2+ exchangers (NCX), JSY N is the flux through the post-synaptic
4 Results
density (PSD), and JV DCC is the flux from voltage dependent calcium channels (VDCCs). Details of how each flux
In this section, we demonstrate the effectiveness of PINNs
term is calculated are given in Appendix A.
in modeling calcium dynamics. All constants used are
given in Table 1 and Table 2. We consider the computational domain (x, t) : [−1, 1] × [0, 1], and use a numeri3 Physics-Informed Neural Networks
cally computed solution as the exact solution (see (3) for
To simulate calcium dynamcis, we incorporate the above details on numerical methods). This is effectively a short
diffusion model into a physics-informed neural network time simulation that is not difficult to numerically simulate,(PINN) (11; 9; 10). Consider a general partial differential but a learned model may be used for long time simulations (100 to 1000 times longer). The network consists ofequation of the form
4 hidden layers of 50 computational neurons each. The
ut + N [u] = 0,(8) network is initialized with a uniform Xavier scheme, and2

Coupled Calcium PINNan Adam optimizer with learning rate of 0.001 is used tominimize the M SE. N = 300 collocation points are randomly generated on both the interior of the domain andthe boundary (50 on each spatial boundary and 50 on theinitial condition). These can be seen in Fig. 1.

tion moves forward in time and the concentration values
become consistent with the diffusion model. By the final
time, the numerical L2 error in the calcium has dropped to
5.83 × 10−4 .
For
TMStreatment purposes,
wemay be interested in what
is occurringat the soma
(cell body) of
the neuron.
Fig 3 shows
Figure 3: Concentrations at soma
the concentration over timeat the soma of our cable neuron. We see that the network
very accurately predicts the concentrations of both calcium and calbindin over time, producing a solution nearly
indistinguishable from the exact solution.

For the biological model, we introduce a constant stimulus
at the soma for the entire simulation time, and a constantvoltage is applied. As can be seen in Appendix A, the gating function for VDCCs requires the solution of two ODEs.
For the purposes of simulating TMS treatment protocols,predetermined voltage data may be provided to review its
effects. Since the gating function only requires voltage andtemporal values, these gating probabilities may be precomputed prior to training the network. We therefore avoid therequirement of modeling additional ODEs, and provide
these gating probabilities explicitly to the network.

5

Conclusions

In this study, we presented a physics-informed neural network (PINN) for modeling a coupled calcium model ina cable neuron. This is significant for the study of repetitive transcranial magnetic stimulation (rTMS), as it allows for long-time simulations with near instantaneous
results. While traditional numerical methods provide accurate solutions, the requirement of small timesteps resultsin significant computation time and hampering the use
for simulating a full rTMS treatment that could last several minutes with voltage pulses in millisecond bursts (13).
Training a PINN on a short timespan allows for feeding
data into the fully trained network and producing a longsimulation.

Figure 1: Domain collocation points
We first lookat the valuesof the initialcalcium andcalbindin concentrations.
In the numerical simulation, the initial concenFigure 2: Initial concentrations
trations are
given as constant values. These may not represent the true steadystate, and they will not satisfy the diffusion model, whichleads to error in the PINN if they are used. We therefore
consider the initial condition as the computed values after
a small number of numerical timesteps (∼5% of the total
timesteps), as the concentrations are approaching values
that satisfy the diffusion model. As the number of numerical timesteps is selected a priori, the values are still notguaranteed to be consistent with the diffusion model, but
this significantly improves the performance of the PINN.

Our results show that PINNs are capable of accurately
capturing the calcium concentration over a cable neuron,
which can be extended over time. While this work is limited to a single cable neuron, using this to create branches,
and therefore a full neuron model, is an important direction
for future work.

Broader Impact
Transcranial Magnetic Stimulation is an important noninvasive treatment for a wide array of neurological disorders.
While the neuronal response to TMS is still poorly understood, simulation tools are providing a way of studying
these effects and how they will benefit patients. A limiting
factor in this research has been the difficulty in developing long-time simulations using traditional numericalmethods. Neural networks provide the ability to provide
long-time simulations with ease that may be applied inboth a research and clinical setting, benefitting researchersand patients alike. Our work is an important step in that
direction, and has the potential to fill this need.

In Fig. 2, we show the results of the network prediction of
the initial concentrations. While there is a small amount oferror, we can see that the network fairly accurately predicts
these values, with a numerical L2 error in the calcium of
1.34 × 10−3 . In Fig. 4 and Fig. 5, we see that this error
is significantly reduced to nearly vanishing as the simula3

Coupled Calcium PINN

Table 1: Values of constants in calcium model and fluxes
Constant

Table 2: Values of constants in calcium model VDCC flux

Value

D
btot
kb−
kb+

Constant
JV DCC

1
5
6
3

R
F
T
z
τk,0
τl,0
Kk
Kl
zk
zl
γk
γl
V1/2,kV1/2,l

JP
ρP
IP
KP

2
1E-2
3
JN

ρN
IN
KN

2E5
1E-4
180
JSY N

jrls

Value

1

Figure 4: Calcium concentrations

Figure 5: Calbindin concentrations

4

8.314
96485
310
2
1.7E-370E-3
1.7E-3
70E-3
2
1
0
0
-21E-3-40E-3

Coupled Calcium PINN

Acknowledgments[10], Physics informed deep learning (part ii):
Data-driven discovery of nonlinear partial differential equations, arXiv preprint arXiv:1711.10566,(2017).
[11] M. R AISSI , P. P ERDIKARIS , AND G. E. K AR NIADAKIS , Physics-informed neural networks: A
deep learning framework for solving forward andinverse problems involving nonlinear partial differential equations, Journal of Computational Physics,
378 (2019), pp. 686–707.
[12] M. H. S AADAT, B. G JORGIEV, L. DAS , AND
G. S ANSAVINI, Neural tangent kernel analysis ofpinn for advection-diffusion equation, arXiv preprint
arXiv:2211.11716, (2022).
[13] S. S HIRINPOUR , N. H ANANEIA , J. ROSADO ,H. T RAN , C. G ALANIS , A. V LACHOS ,P. J EDLICKA , G. Q UEISSER , AND A. O PITZ,
Multi-scale modeling toolbox for single neuron
and subcellular activity under transcranial magnetic stimulation, Brain Stimulation, 14 (2021),
pp. 1470–1482.
[14] V. S ITZMANN , J. M ARTEL , A. B ERGMAN , D. L IN DELL , AND G. W ETZSTEIN , Implicit neural representations with periodic activation functions, inAdvances in Neural Information Processing Systems,
vol. 33, 2020, pp. 7462–7473.
[15] W. S NYDER , I. T EZUAR , AND C. R. W ENTLAND,
Domain decomposition-based coupling of physicsinformed neural networks via the schwarz alternating
method, arXiv preprint arXiv:2311.00224, (2023).

This work was supported by NIH grant R01EB034143.

References[1] S. A BRAHAM , J. K INNISON , Z. M. M IKSIS ,D. P OSTER , S. YOU , J. D. H AUENSTEIN , AND
W. S CHEIRER, Efficient hyperparameter optimization for atr using homotopy parametrization, in Automatic Target Recognition XXXIII, R. I. Hammoud,T. L. Overman, and A. Mahalanobis, eds., vol. 12521,
2023, p. 1252107.
[2] A. T. BARKER , R. JALINOUS , AND I. L. F REESTON,
Non-invasive magnetic stimulation of human motor
cortex, The Lancet, 325 (1985), pp. 1106–1007.
[3] P. B OROLE, CalciumSim: Simulator for calcium dynamics on neuron graphs using dimensionally reduced model, 2022.
[4] P. R. B OROLE , J. M. ROSADO , M. N EAL , AND
G. Q UEISSER, Neuronal resilience and calcium signaling pathways in the context of synapse loss andcalcium leaks: A computational modeling study andimplications for alzheimer’s disease, SIAM Journal
on Applied Mathematics, 83 (2023), pp. 2418–2442.
[5] S. G REIN , M. S TEPNIEWSKI , S. R EITER , M. M.
K NODEL , AND G. Q UEISSER, 1d-3d hybrid
modeling-from multi-compartment models to full resolution models in space and time, Frontiers in Neuroinformatics, 8 (2014).
[6] Q. G UAN AND G. Q UEISSER, Modeling calcium dynamics in neurons with endoplasmic reticulum: Existence, uniqueness and an implicit–explicit finite element scheme, Communications in Nonlinear Scienceand Numerical Simulation, 109 (2022), p. 106354.
[7] M. H ALLETT, Transcranial magnetic stimulation: A
primer, Neuron, 55 (2007), pp. 187–199.
[8] J.-P. L EFAUCHEUR , N. A NDRÉ -O BADIA , A. A N TAL , S. S. AYACHE , C. BAEKEN , D. H. B EN NINGER , R. M. C ANTELLO , M. C INCOTTA , M. DE
C ARVALHO , D. D E R IDDER , H. D EVANNE , V. D I
L AZZARO , S. R. F ILIPOVI Ć , F. C. H UMMEL , S. K.
JÄÄSKELÄINEN , V. K. K IMISKIDIS , G. KOCH ,B. L ANGGUTH , T. N YFFELER , A. O LIVIERO ,F. PADBERG , E. P OULET, S. ROSSI , P. M. ROSSINI ,J. C. ROTHWELL , C. S CHÖNFELDT-L ECUONA ,H. R. S IEBNER , C. W. S LOTEMA , C. J. S TAGG ,J. VALLS -S OLE , U. Z IEMANN , W. PAULUS , AND
L. G ARCIA -L ARREA, Evidence-based guidelines on
the therapeutic use of repetitive transcranial magnetic stimulation (rtms), Clinical Neurophysiology,
125 (2014), pp. 2150–2206.
[9] M. R AISSI , P. P ERDIKARIS , AND G. E. K ARNI ADAKIS , Physics informed deep learning (part i):
Data-driven solutions of nonlinear partial differential
equations, arXiv preprint arXiv:1711.10561, (2017).
5

Coupled Calcium PINN

Appendix A: Membrane Flux Components

where G(V, t) ∈ [0, 1] is the gating function and
F (V, ∆[Ca2+ ]) is the flux function. The difference beThe plasma membrane flux in this study can be separated tween cytoplasmic and extracellular ion concentration is
given byas
JP M = −JP − JN + JSY N + JV DCC .

∆[Ca2+ ] = cc − co ,(14)

In this section, we give a brief description of each compo- and the flux function is given by the Goldman-HodgkinKatz equation,nent as described in (4).
PMCA Pumps A second-order Hill equation is used
2 2
to model the plasma membrane Ca2+ current of PMCA F (V, ∆[Ca2+ ]) = pCa2+ V z F · cc − co exp(−zF V /RT )RT
1 − exp(zF V /RT )pumps,(20)
IP c2c,KP2 + c2cwhere R is the universal gas constant, F is Faraday’s con(15) stant, T is temperature in Kelvin, p 2+ is the permeability
Caof Ca2+ ions through the channels, and z is the valence
where ρP is the density of PMCA pumps on the plasma of a Ca2+ ion. The gating function for N-type VDCCs ismembrane, IP is the single channel Ca2+ current, cc is the given bycytosolic Ca2+ concentration, and KP is the measure of
Ca2+ affinity.
G(V, t) = k(V, t)l2 (V, t),(21)
JP (cc ) = ρP ·NCX Exchangers A second-order Hill equation is used where the gating functions k(·) and l(·) satisfy the two
to model the plasma membrane Ca2+ current of NCX ex- ODEschangers,∂kk∞ − kl∞ − l
∂l
=
=
and
,
IN cc
∂t
τk
∂t
τl
JN (cc ) = ρN ·(16)
KN + cc
where
where ρN is the density of NCX exchangers on the plasma
membrane, IN is the single channel Ca2+ current, cc is theαk (V )
αl (V )
cytosolic Ca2+ concentration, and KN is the measure of
k∞ =
;
l∞ =
,
2+
α
(V
)
+
β
(V
)
α
(V
) + βl (V )kkl
Ca affinity.
Synaptic Influx Calcium influx is modeled by a linearly anddecreasing function:t − t0JSY N = jrls · 1 −
λt0 (t)
τrls

τk =
(17)

1
+ τk,0 ;αk + βkzk γk (V − V1/2,k )FRT−zk (1 − γk )(V − V1/2,k )F
βk (V ) = Kk exp
RTzl γl (V − V1/2,l )F
αl (V ) = Kl expRT−zl (1 − γl )(V − V1/2,l )F
βl (V ) = Kl expRT

λt0 (t) =1,0,t ∈ [t0 , t0 + τrls ]otherwise.

αk (V ) = Kk exp(18)

VDCCs For VDCCs, we use a Borg-Graham model. The
Ca2+ current is given by
JV DCC (V, cc , t) = G(V, t)F (V, ∆[Ca2+ ]),1+ τl,0 .
αl + βlThe rate functions are defined aswhere jrls is the current density and τrls is the associatedtime constant. The function λt0 (t) is given by:


τl =

(19)

6