Modeling excitable cells with the EMI equations:

spectral analysis and iterative solution strategy

arXiv:2308.12145v1 [math.NA] 23 Aug 2023

Pietro Benedusi, Paola Ferrari, Marie E. Rognes,

Stefano Serra-Capizzano

Abstract
In this work, we are interested in solving large linear systems stem-
ming from the Extra-Membrane-Intra (EMI) model, which is employed for
simulating excitable tissues at a cellular scale. After setting the related
systems of partial differential equations (PDEs) equipped with proper
boundary conditions, we provide numerical approximation schemes for
the EMI PDEs and focus on the resulting large linear systems. We first
give a relatively complete spectral analysis using tools from the theory
of Generalized Locally Toeplitz matrix sequences. The obtained spec-
tral information is used for designing appropriate (preconditioned) Krylov
solvers. We show, through numerical experiments, that the presented so-
lution strategy is robust w.r.t. problem and discretization parameters,
efficient and scalable.

1 Introduction
The EMI (Extra, Intra, Membrane) model, also known as the cell-by-cell model,
is a numerical building block in computational electrophysiology, employed to
simulate excitable tissues at a cellular scale. With respect to homogenized
models, such as the well-known mono and bidomain equations, the EMI model
explicitly resolves cells morphologies, enabling detailed biological simulations.
For example, inhomogeneities in ionic channels along the membrane, as observed
in myelinated neuronal axons [12], can be described within the EMI model.
The main areas of application for the EMI model are computational cardiol-
ogy and neuroscience, where spreading excitation by electrical signalling plays
a crucial role [13, 17, 25, 18, 24, 14, 33, 29]. We refer to [34] for an exhaustive
description of the EMI model, its derivation, and relevant applications.
This work considers a Galerkin-type approximation of the underlying system
of partial differential equations (PDEs). A finite element discretization leads to
block-structured linear systems of possibly large dimensions. The use of ad hoc
tools developed in the theory of Generalized Locally Toeplitz matrix sequences
allows us to give a quite complete picture of the spectral features (the eigenvalue
distribution) of the resulting matrices and matrix sequences.

1
 The latter information is crucial to understand the algebraic properties of
such discrete objects and to design iterative solution strategies, enabling large-
scale models, as in [10, 9]. More precisely, spectral analysis is employed for
designing tailored preconditioning strategies for Krylov methods. The resulting
efficient and robust solution strategies are theoretically studied and numerically
tested.
The current paper is organized as follows. In Section 2 the continuous prob-
lem is introduced together with possible generalizations. Section 3 considers a
basic Galerkin strategy for approximating the examined problem. The spectral
analysis is given in Section 4 in terms of distribution results and degenerating
eigenspaces. Specifically, Subsection 4.1 lays out the foundational theories and
concepts necessary for understanding the distribution of the entire system (pre-
sented in a general form in Subsection 4.2) and the specific stiffness matrices
and matrix sequences (discussed in Subsection 4.3). These findings are then em-
ployed in Subsection 4.4 for proposing a preconditioning strategy. Moreover, in
Section 5 we report and critically discuss the numerical experiments and Section
6 contains conclusions and a list of relevant open problems.

2 The EMI problem

We introduce the partial differential equations characterizing the EMI model.
When comparing to the homogenized models, we observe that the novelty of
the EMI approach consists in the fact that the cellular membrane Γ is explicitly
represented, as well as the intra- and extra-cellular quantities, denoted with
subscript i and e respectively.
More in detail, given a domain Ω = Ωi ∪ Ωe ∪ Γ ⊂ Rd , typically with
d ∈ {2, 3}, ∂Ω = ∂Ωe /Γ, and ∂Ωi = Γ, we consider the following time-discrete
problem for the intra- and extra-cellular potentials ui , ue , and for the membrane
current Im :

− ∇ · (σe ∇ue (x)) = 0 for x ∈ Ωe , (1)

− ∇ · (σi ∇ui (x)) = 0 for x ∈ Ωi , (2)
σe ∇ue (x) · ne = −σi ∇ui (x) · ni ≡ Im (x) for x ∈ Γ, (3)
ui (x) − ue (x) = v(x) for x ∈ Γ, (4)
v(x) − τ Im (x) = f (x) for x ∈ Γ, (5)

with σe , σi , τ ∈ R+ , ni (resp. ne ) is the outer normal on ∂Ωi (resp. ∂Ωe )

and f ∈ L2 (Γ) is known. We can close the EMI problem with homogeneous
boundary conditions:

ue (x) = 0 for x ∈ ∂ΩD , (6)

σe ∇ue (x) · ne = 0 for x ∈ ∂ΩN , (7)

with ∂Ω = ∂ΩD ∪ ∂ΩN . In the case of a pure Neumann problem, with ∂Ω =

∂ΩN , uniqueness must be enforced through an additional constraint (e.g. a

2
condition on the integral of ue ). We are imposing homogeneous boundary con-
ditions for simplicity; the inhomogeneous case reduces to the homogeneous case
by considering the lifting of the boundary data.
Essentially, the EMI problem consists of two homogeneous Poisson problems
coupled at the interface Γ through a Robin-type condition (3)-(5), depending on
the current Im . The EMI problem can be considered a mixed-dimensional prob-
lem since the unknowns of interest are defined on sets of different dimensionality,
d for ui and ue and d − 1 for Im .
It is worth noticing that possible dynamics originate in the membrane via
the source term f . In particular, eq. (5) is obtained by discretizing point-wise
the capacitor current-voltage relation, a time-dependent ODE with t ∈ (0, T ]
given a final time T > 0:

∂v(t) −1
= Cm (Im (t) − Iion (v, t)),
∂t
v(0) = v0 ,

where the positive constant Cm is a capacitance and the ionic current Iion is a
reaction term. An implicit (resp. explicit) integration of Im (resp. Iion ), with
time step ∆t > 0, results in eq. (5) with

f (x) = v0 (x) − τ Iion (v0 (x)), (8)

−1
and τ = Cm ∆t.
The model in eq. (1)-(5) can describe a single cell or multiple disjoint cells,
SNcell j
i.e. with Ωi = j=1 Ωi , in an extracellular media, cf. Figure 1. Extending
model (1)-(5) to multiple cells, with possibly common membranes (a.k.a. gap
junctions), is straightforward [23].

3 Weak formulation and discrete operators

The EMI problem can be weakly formulated in various ways, depending on
the unknowns of interest. We refer to [34] for a broad discussion on various
formulations (including the so-called mixed ones). As it could be expected
from the structure of (1)-(5), all formulations and corresponding discretizations
give rise to block operators, with different blocks corresponding to Ωi , Ωe , and
possibly Γ.
We use the so-called single-dimensional formulation and the corresponding
discrete operators. In this setting, the weak form depends only on bulk quanti-
ties ui and ue since the current term Im is replaced by:

Im (x) = τ −1 (ui (x) − ue (x) − f (x)),

according to equations (4)-(5). Let us remark that “single” refers to the previous
substitution, eliminating the variable defined in Γ; the overall EMI problem is
still in multiple dimensions, in the sense that d ≥ 1.

3
 After substituting the expression for Im in (3), assuming the solution ur , for
r ∈ {i, e}, to be sufficiently regular over Ωr , we multiply the PDEs in (1)-(2) by
test functions vr (x) ∈ Vr (Ωr ), with Vr a sufficiently regular Hilbert space with
elements satisfying the boundary conditions in (6)-(7); in practice H 1 (Ωi ) and
H 1 (Ωe ), would be a standard choice. After integrating over Ωr and applying
integration by parts, using the normal flux definition (3), the weak EMI problem
reads: find ui ∈ Vi (Ωi ) and ue ∈ Ve (Ωe ) such that

Z Z Z Z
τ σe ∇ue · ∇ve dx + ue ve ds −
ui ve ds = − f ve ds, (9)
Ωe Γ Γ Γ
Z Z Z Z
τ σi ∇ui · ∇vi dx + ui vi ds − ue vi ds = f vi ds, (10)
Ωi Γ Γ Γ

for all test functions ve ∈ Ve (Ωe ) and vi ∈ Vi (Ωi ). We refer to [34, Section 6.2.1]
for boundedness and coercivity results for this formulation.
For each subdomain Ωr , with r ∈ {i, e}, we construct a conforming tassel-
lation Tr . We then introduce a yet unspecified discretization via finite element
basis functions (e.g. Lagrangian elements of order p ∈ N on a regular grid) for
Ve and Vi :
   
Ve,h = span {ϕej }N e
j=1 , Vi,h = span {ϕ i Ni
}
j j=1 ,

with Ne , Ni ∈ N denoting the number of degrees of freedom in the correspond-

ing subdomains. We can further decompose Nr = Nr,in + Nr,Γ , i.e. further
differentiating between internal and membrane degrees of freedom, with Nr,Γ
basis functions ϕrj having support intersecting Γ, cf. Figure 1. In the numerical
experiments, we will consider matching Tr on the interface Γ and the same p for
Ve,h and Vi,h ; nevertheless Section 3 and Section 4 are developed in a general
setting, so that the theory is ready also for potential extensions.
From (9)-(10) we define the following discrete operators: intra- and extra-
Laplacians
Z Ne
Ae = σe ∇ϕej (x) · ∇ϕek (x) dx ∈ RNe ×Ne , (11)
Ωe j,k=1
Z Ni
Ai = σi ∇ϕij (x) · ∇ϕik (x) dx ∈ RNi ×Ni , (12)
Ωi j,k=1

membrane mass matrices:

Z Ne
Me = ϕej (x)ϕek (x) ds ∈ RNe ×Ne , (13)
Γ j,k=1
Z Ni
Mi = ϕij (x)ϕik (x) ds ∈ RNi ×Ni , (14)
Γ j,k=1

4
 ∂Ω e,in

Ωe ∂Ω Vi,h (Ωi ) e, Γ
Ω1i
i,in
Γ
Ω2i
Γ i, Γ
Ve,h (Ωe )

Figure 1: Left: example of EMI geometry with Ωi = Ω1i ∪ Ω2i . Right: example
of the discretization setting for d = 2 on a regular grid. We notice that interior
and exterior nodes overlap in Γ. For p = 1, in this case, we have Ni = 9, Ne = 24
with Ni,in = 1, Ne,in = 16, and NΓ = 8. Corresponding labels are reported.

and the coupling matrix

Z (Ne ,Ni )
Tei = ϕij (x)ϕek (x) ds ∈ RNe ×Ni . (15)
Γ (j,k)=(1,1)

Finally, we write the linear system of size n = Ne + Ni corresponding to (9)-(10)

as     
τ A e + Me Tei ue f
T = e ⇐⇒ An u = f , (16)
Tei τ Ai + Mi ui fi
with ur , fr ∈ RNr the unknowns and the right hand side corresponding to Ωr
and r ∈ {i, e}. We remark that the operator An is symmetric and positive
definite1 . Making the degrees of freedom in Γ denoted explicitly (with reference
to Figure 1), i.e. ui,Γ , ue,Γ ∈ RNΓ , as well as the interior ones, i.e. ui,in ∈
RNi −NΓ and ue,in ∈ RNe −NΓ , we can rewrite more extensively system (16) as

τ Ain τ Ain,Γ
    
e e 0 0 ue,in 0
τ AΓ,in
e τ AΓe + MeΓ 0 Γ
Tei   ue,Γ  fΓ,e 
 in,Γ  
 = , (17)
 0 0 τ Aini τ Ai ui,in   0 
0 Γ T
(Tei ) τ AΓ,in
i τ AΓi + MiΓ ui,Γ fΓ,i

with
τ Ain τ Ain,Γ
 
e e 0 0
τ AΓ,in
e τ AΓe + MeΓ 0 Γ
Tei 
An = 
 in,Γ  .
 (18)
0 0 τ Aini τ Ai
0 Γ T
(Tei ) τ AΓ,in
i τ AΓi + MiΓ
1 Dirichlet boundary conditions can be imposed to the right-hand side to enforce symmetry.

5
For forthcoming use, we define the bulk matrix
 in 
τ Ae 0 0 0
 0 0 0 0
Ãn =  ,
 0 0 τ Ain
i 0
0 0 0 0

where all membrane terms are zeroed.

4 Spectral analysis
In this section, we study the spectral distribution of the matrix sequence {An }n
under various assumptions for determining the global behaviour of the eigen-
values of An as the matrix size n tends to infinity. The spectral distribution is
given by a smooth function called the (spectral) symbol as it is customary in
the Toeplitz and Generalized Locally Toeplitz (GLT) setting [20, 21, 8, 7].
First, we give the formal definition of Toeplitz structures, eigenvalue (spec-
tral) and singular value distribution, the basic tools that we use from the GLT
theory, and finally, we provide the specific analysis of our matrix sequences
under a variety of assumptions.

4.1 Toeplitz structures, spectral symbol, GLT tools

We initially formalize the definition of block Toeplitz and circulant sequences
associated with a matrix-valued Lebesgue integrable function. Then, we pro-
vide the notion of eigenvalue (spectral) and singular value distribution, and we
introduce the basic tools taken from the GLT theory.

Definition 4.1 [Toeplitz, block-Toeplitz, multilevel Toeplitz matrices] A finite-

dimensional or infinite-dimensional Toeplitz matrix is a matrix that has constant
elements along each descending diagonal from left to right, namely,

a−1 a−2 · · · a−n+1

 
  a0
a0 a−1 a−2 · · ·  .. 
 a1 a a · · ·  a1 a 0 a −1 . 
0 −1   
T =

. .

, T =
 . . . . .
. 
(19)
 a2 a .. . .  n  a2
  a1 . . . 
1 
   .
.. .. .. ..  ..

. . . . a1 a0 a−1 
an−1 · · · a2 a1 a0

As the first equation in (19) indicates, the matrix T can be infinite-dimensional.

A finite-dimensional Toeplitz matrix of dimension n is denoted as Tn . We also
consider sequences of Toeplitz matrices as a function of their dimension, denoted
as {Tn }n , n = 1, 2, · · · , ∞.
In general, the entries ak , k ∈ Z, can be matrices (blocks) themselves, defin-
ing Tn as a block Toeplitz matrix. Thus, in the block case, ak , k ∈ {−n +
1, · · · , n−1} are blocks of size s1 ×s2 , the subindex of Tn is the number of blocks

6
in the Toeplitz matrix while N1 × N2 is the size of the matrix with N1 = n s1 ,
N2 = n s2 . To be specific, we use also the notation XN1 ,N2 = Tn . A special case
of block-Toeplitz matrices is the class of two- and multilevel block Toeplitz ma-
trices, where the blocks are Toeplitz (or multilevel Toeplitz) matrices themselves.
The standard Toeplitz matrices are sometimes addressed as unilevel Toeplitz.

Definition 4.2 [Toeplitz sequences (generating function of)] Denote by f a

d-variate complex-valued integrable function, defined over the domain Qd =
[−π, π]d , d ≥ 1, with d-dimensional Lebesgue measure µd (Qd ) = (2π)d . Denote
by fk the Fourier coefficients of f ,
Z
1
fk = f (θ)e−i (k,θ) dθ, k = (k1 , · · · , kd ) ∈ Zd , i2 = −1,
(2π)d
Qd

Pd
where θ = (θ1 , · · · , θd ), (k, θ) = j=1 kj θj , n = (n1 , · · · , nd ), and N (n) =
n1 × · · · × nd . By following the multi-index notation in [35][Section 6], with each
f we can associate a sequence of Toeplitz matrices {Tn }n , where

Tn = {fk−ℓ }nk,ℓ=eT ∈ CN (n)×N (n) ,

e = [1, 1, · · · , 1] ∈ Nd . For d = 1
 
f0 f−1 ··· f−n+1
 .. .. 
 f1 f0 . . 
Tn = 
 .
,
 .. .. 
. f0 f−1 
fn−1 ··· f1 f0

or for d = 2, i.e. the two-level case, and for example n = (2, 3), we have
 
  fk,0 fk,−1 fk,−2
F0 F−1
Tn = , Fk = fk,1 fk,0 fk,−1  , k = 0, ±1.
F1 F0
fk,2 fk,1 fk,0

The function f is referred to as the generating function (or the symbol of ) Tn .

Using a more compact notation, we say that the function f is the generating
function of the whole sequence {Tn }n and we write Tn = Tn (f ).
If f is d-variate, Cs1 ×s2 matrix-valued, and integrable over Qd , d, s1 , s2 ≥ 1,
i.e. f ∈ L1 (Q2 , s1 × s2 ), then we can define the Fourier coefficients of f in
the same way (now fk is a matrix of size s1 × s2 ) and consequently Tn =
{fk−ℓ }nk,ℓ=eT ∈ Cs1 N (n)×s2 N (n) , then Tn is a d-level block Toeplitz matrix ac-
cording to Definition 4.1. If s1 = s2 = s then we write f ∈ L1 (Q2 , s).
As in the scalar case, the function f is referred to as the generating func-
tion of Tn . We say that the function f is the generating function of the whole
sequence {Tn }n , and we use the notation Tn = Tn (f ).

7
Definition 4.3 Let f : D → Cs×s be a measurable matrix-valued function with
eigenvalues λi (f ) and singular values σi (f ), i = 1, . . . , s. Assume that D ⊂ Rd
is Lebesgue measurable with positive and finite Lebesgue measure µd (D). Assume
that {An }n is a sequence of matrices such that dim(An ) = dn → ∞, as n → ∞
and with eigenvalues λj (An ) and singular values σj (An ), j = 1, . . . , dn .

• We say that {An }n is distributed as f over D in the sense of the eigen-

values, and we write {An }n ∼λ (f, D), if
dn Z s
1 X 1 1X
lim F (λj (An )) = F (λi (f (t))) dt, (20)
n→∞ dn µd (D) D s i=1
j=1

for every continuous function F with compact support. In this case, we

say that f is the spectral symbol of {An }n .
• We say that {An }n is distributed as f over D in the sense of the singular
values, and we write {An }n ∼σ (f, D), if
dn Z s
1 X 1 1X
lim F (σj (An )) = F (σi (f (t))) dt, (21)
n→∞ dn µd (D) D s i=1
j=1

for every continuous function F with compact support. In this case, we

say that f is the singular value symbol of {An }n .
• The notion {An }n ∼σ (f, D) applies also in the rectangular case where
f is Cs1 ×s2 matrix-valued. In such a case the parameter s in formula
(21) has to be replaced by the minimum between s1 and s2 : furthermore
(1) (2) (1)
An ∈ Cdn ×dn with dn in formula (21) being the minimum between dn
(2)
and dn . Of course the notion of eigenvalue distribution does not apply in
a rectangular setting.

Throughout the paper, when the domain can be easily inferred from the context,
we replace the notation {An }n ∼λ,σ (f, D) with {An }n ∼λ,σ f .

Remark 4.1 If f is smooth enough, an informal interpretation of the limit

relation (20) (resp. (21)) is that when n is sufficiently large, the eigenvalues
(resp. singular values) of An can be subdivided into s different subsets of the
same cardinality. Then dn /s eigenvalues (resp. singular values) of An can be
approximated by a sampling of λ1 (f ) (resp. σ1 (f )) on a uniform equispaced
grid of the domain D, and so on until the last dn /s eigenvalues (resp. singular
values), which can be approximated by an equispaced sampling of λs (f ) (resp.
σs (f )) in the domain D.

Remark 4.2 We say that An is zero-distributed if An ∼λ 0. Of course, if the

eigenvalues of An tend all to zero for n → ∞, then this is sufficient to claim
that An ∼λ 0.

8
For Toeplitz matrix sequences, the following theorem due to Tilli holds, which
generalizes previous research along the last 100 years by Szegő, Widom, Avram,
Parter, Tyrtyshnikov, and Zamarashkin (see [21, 7] and references therein).
Theorem 4.1 [32] Let f ∈ L1 (Qd , s1 × s2 ), then {Tn (f )}n ∼σ (f, Qd ). If
s1 = s2 = s and if f is a Hermitian matrix-valued function, then {Tn (f )}n ∼λ
(f, Qd ).
The following theorem is useful for computing the spectral distribution of
a sequence of Hermitian matrices. For the related proof, see [27, Theorem 4.3]
and [28, Theorem 8]. Here, the conjugate transpose of the matrix X is denoted
by X ∗ .
Theorem 4.2 [27, Theorem 4.3] Let {An }n be a sequence of matrices, with
An Hermitian of size dn , and let {Pn }n be a sequence such that Pn ∈ Cdn ×δn ,
Pn∗ Pn = Iδn , δn ≤ dn and δn /dn → 1 as n → ∞. Then {An }n ∼λ f if and only
if {Pn∗ An Pn }n ∼λ f .
With the notations of the result above, the matrix sequence {Pn∗ An Pn }n
is called a compression of {An }n and the single matrix Pn∗ An Pn is called a
compression of An .
In what follows we take into account a crucial fact that often is neglected: the
generating function of a Toeplitz matrix sequence and even more the spectral
symbol of a given matrix sequence is not unique, except for the trivial case of
either a constant generating function or a constant spectral symbol. In fact,
here we report and generalize Remark 1.3 at p. 76 in [15] and the discussion
below Theorem 3 at p. 8 in [16].
Remark 4.3 [Multiple Toeplitz generating functions and multiple spectral sym-
bols] Let n be a multiple of k and consider the Toeplitz matrix Xn = Tn (f ). We
can view Xn as a block-Toeplitz with blocks of size k ×k such that Xn = Tn (f ) =
T nk (f [k] ). In our specific context, we have

f [k] (θ) = Tk (f ) − e1 eTk ei θ

− ek eT1 e−i θ

with ej , j = 1, . . . , k, being the canonical basis of Ck .

As an example, let n be even and consider the function f (θ) = 2 − 2 cos(θ).
According to Definition 4.2, Xn = Tn (f ) but we can also view the matrix as
Xn = T n2 (f [2] ), namely
   
2 −1 0 A0 A−1 0
−1 . . . . . .  A1 . . . . . .
   
 = T n (f [2] ),

Xn = Tn (f ) =  
=
 
. .. . . . −1  . .. . .. A 
 2

−1

0 −1 2 0 A1 A0
where    
2 −1 0 −1
A0 = , A1 = AT−1 = .
−1 2 0 0

9
Thus, f [2] = A0 + A1 eiθ + AT1 e−iθ .

It is clear that analogous multiple representations of Toeplitz matrices hold

for every function f . As a consequence, taking into consideration Theorem 4.1,
we have simultaneously {Tn (f )}n ∼λ (f, Qd ) and {Tn (f )}n ∼λ (f [k] , Qd ) for
any fixed k independent of n. It should be also observed that a situation in
which k does not divide n is not an issue thanks to the compression argument
in Theorem 4.2.
More generally, we can safely claim that the spectral symbol in the Weyl
sense of Definition 4.3 is far from unique and in fact any rearrangement is still
a symbol (see [16, 6]). A simple case is given by standard Toeplitz sequences
{Tn (f )}n , with f real-valued and even that is f (θ) = f (−θ) almost everywhere,
θ ∈ Q. In that case
Z π
1 π
Z
1
F (f (θ)) dθ = F (f (θ)) dθ, (22)
2π −π π 0
due to the even character of f , and hence it is also true that {Tn (f )}n ∼λ
(f, Q+ ), Q+ = (0, π). A general analysis of these concepts via rearrangement
theory can be found in [6] and references therein.

4.2 Symbol analysis

In this section, we state and prove three results which hold under reasonable
assumptions. The first has the maximal generality, while the second and the
third can be viewed as special cases of the first. Let us remind that Ni , Ne and
NΓ denote the number of intra-, extra- and membrane degrees of freedom.
Theorem 4.3 Assume that
NΓ = o(min{Ni , Ne }) for NΓ , Ni , Ne → ∞,
Ni
lim = r ∈ (0, 1).
Ni ,Ne →∞ Ne + Ni

Assume that
{τ Ain e e
e }n ∼λ (f , D ), (23)
{τ Ain
i }n
i
∼λ (f , D ),i
(24)

De ⊂ Rke , Di ⊂ Rki , given f e and f i extra- and intra- spectral symbols. It

follows that
{Ãn }n , {An }n ∼λ (g, [0, 1] × Di × De )
where g(x, ti , te ) = f i (ti )ψ[0,r] (x) + f e (te )ψ(r,1] (x), x ∈ [0, 1], ti ∈ Di , te ∈ De ,
and
τ Ain,Γ
 
0 e 0 0
τ AeΓ,in Γ Γ Γ
τ A e + Me 0 Tei 
 ∼λ 0,
τ Ain,Γ

 0 0 0 i

Γ T Γ,in Γ Γ
0 (Tei ) τ Ai τ A i + Mi
ψZ denoting the characteristic function of the set Z.

10
Proof
First, we observe that because of the Galerkin approach and the structure
of the equations, all the involved matrices are real and symmetric. Hence the
symbols f e , f i are necessarily Hermitian matrix-valued. Without loss of gen-
erality we assume that both f e , f i take values into Cs×s (the case where f e ,
f i have different matrix sizes is treated in Remark 4.4, for the sake of nota-
tional simplicity). Because of equations (23) and (24), setting Ne′ = Ne − NΓ ,
Ni′ = Ni − NΓ , and taking into account Definition 4.3, we have
′
Ne Z s
1 X 1 1 X
lim F (λj (τ Ain
e )) = F (λm (f e (te ))) dte , (25)
Ne →∞ Ne′ µke (De ) De s m=1
j=1
′
Ni Z s
1 X 1 1 X
lim F (λj (τ Ain
i )) = F (λm (f i (ti ))) dti , (26)
Ni →∞ Ni′ µki (Di ) Di s m=1
j=1

for any continuous function F with bounded support. Now we build the 2 × 2
block matrix
X[Ne , Ni , NΓ ] = diag(τ Ain in
e , τ Ai )

and we consider again a generic continuous function F with bounded support.

By the block diagonal structure of X[Ne , Ni , NΓ ], we infer

′ ′
Ne +Ni
1 X
∆(Ne , Ni , NΓ , F ) = ′ ′ F (λj (X[Ne , Ni , NΓ ]))
Ne + Ni j=1
′ ′
Ne Ni
X F (λj (τ Ain )) e
X F (λj (τ Ain )) i
= +
j=1
Ne′ + Ni′ j=1
Ne′ + Ni′
Ne′ Ni ′
Ne′ X F (λj (τ Ain ))
e Ni′ X F (λj (τ Ain
i ))
= ′ + ′ .
Ne′ + Ni j=1 Ne′ Ne′ + Ni j=1 Ni′

As a consequence, taking the limit as Ni , Ne → ∞, using the assumption

NΓ = o(min{Ni , Ne }),
Ni
lim = r ∈ (0, 1),
Ni ,Ne →∞ Ne + Ni
we obtain that the limit of ∆(Ne , Ni , NΓ , F ) exists and it is equal to
s s
1−r
Z Z
r 1 X 1 X
F (λm (f e (te ))) dte + F (λm (f i (ti ))) dti .
µke (De ) D e s m=1
µki (D i)
D i s m=1
(27)
With regard to Definition 4.3, the quantity (27) does not look like the right-hand
side of (20), since we see two different terms. The difficulty can be overcome by

11
enlarging the space with a fictitious domain [0, 1] and interpreting the sum in
(27) as the global integral involving a step function. In reality, we rewrite (27)
as
Z 1 Z s
1 1 X
dx F (λm (f e (te ))ψ[0,r] (x) + λm (f i (ti ))ψ(r,1] (x)) dte dti ,
µq (D) 0 D e ×D i s m=1

with D = [0, 1] × Di × De , q = ke + ki + 1, µq (D) = µke (De )µki (Di ). Hence the

proof of the relation

{X[Ne , Ni , NΓ ]}n ∼λ (g, [0, 1] × Di × De )

is complete. Now we observe that X[Ne , Ni , NΓ ] is a compression of Ãn (cf.

definition after Theorem 4.2): therefore by exploiting again the hypothesis NΓ =
o(min{Ni , Ne }) which implies NΓ = o(n), by Theorem 4.2, we deduce

{Ãn }n ∼λ (g, [0, 1] × Di × De ).

Finally in Exercise 5.3 in the book [20] it is proven the following: if {Vn }n ∼λ f
(not necessarily Toeplitz or GLT), {An = Vn +Yn }n with {Yn }n zero-distributed
and An , Yn Hermitian for every dimension, then {An }n ∼λ f . This setting is
exactly our setting since by direct inspection the rank of

τ Ain,Γ
 
0 e 0 0
τ AΓ,in
e τ AΓe + MeΓ 0 Γ
Tei 
Rn =  0 in,Γ 

0 0 τ Ai
0 Γ T
(Tei ) τ AΓ,in
i τ AΓi + MiΓ

is bounded by 4NΓ . Hence the number of nonzero eigenvalues of the latter

matrix cannot exceed 4NΓ = o(n), n = Ni + Ne , and therefore the related
matrix sequence is zero-distributed in the eigenvalue sense i.e. {Rn }n ∼λ 0.
Since An = Ãn + Rn , the application of the statement in [20, Exercise 5.3]
implies directly {An }n ∼λ (g, [0, 1] × Di × De ) and the proof is concluded. •

Remark 4.4 The case where f e and f i have different matrix sizes would com-
plicate the derivations in the proof of the above theorem. However, in this case,
the argument is simple and relies on the non-uniqueness of the spectral symbol:
for a discussion on the matter refer to [16, 6] and Remark 4.3.

The following two corollaries simplify the statement of Theorem 4.3, under
special assumptions which are satisfied for a few basic discretization schemes
and when dealing with elementary domains.

Corollary 4.1 Assume that

NΓ = o(min{Ni , Ne }) for NΓ , Ni , Ne → ∞,
Ni
lim = 1/2.
Ni ,Ne →∞ Ne + Ni

12
Assume that
{τ Ain e
e }n ∼λ (f , D), (28)
{τ Ain
i }n
i
∼λ (f , D), (29)
D ⊂ Rk . It follows that
fe
  
0
{Ãn }n , {An }n ∼λ ,D
0 fi
and
τ Ain,Γ
 
0 e 0 0
τ AΓ,in τ AΓe + MeΓ 0 Γ
Tei
 e

in,Γ  ∼λ 0.

 0 0 0 τ Ai
0 Γ T
(Tei ) τ AΓ,in
i τ AΓi + MiΓ
Proof
Notice that the writing
 e 
f 0
{Ãn }n , {An }n ∼λ
0 fi
and
{Ãn }n , {An }n ∼λ f i (θ)ψ[0,c] (x) + f e (θ)ψ(c,1] (x).
are equivalent for c = 1/2.
•
Corollary 4.2 Assume that
NΓ = o(min{Ni , Ne }) for NΓ , Ni , Ne → ∞,
Ni
lim = r ∈ [0, 1].
Ni ,Ne →∞ Ne + Ni

Assume that
{τ Ain
e }n ∼λ (f, D), (30)
{τ Ain
i }n ∼λ (f, D), (31)
D ⊂ Rk . It follows that
{Ãn }n , {An }n ∼λ (f, D)
and
τ Ain,Γ
 
0 e 0 0
τ AeΓ,in τ AΓe + MeΓ 0 Γ
Tei 
in,Γ  ∼λ 0.
 
 0 0 0 τ Ai
0 Γ T
(Tei ) τ AΓ,in
i τ AΓi + MiΓ

Proof The proof of the relation {Ãn }n , {An }n ∼λ (f, D) follows directly from
the limit displayed in (27), after replacing f e , f i with f and necessarily De , Di
with D. The rest is obtained verbatim as in Theorem 4.3.
•

13
4.3 Analysis of stiffness matrix sequences
The spectral analysis of stiffness matrix sequences is crucial to understand the
global spectral behaviour of EMI matrices and, in turn, to justify our precon-
ditioning approach. According to [22, Section 5], given a bi-dimensional square
domain, the sequence of corresponding stiffness matrices {Asquare
n }n obtained
from the linear Q1 Lagrangian finite element method2 (FEM) has the following
spectral distribution:
{Asquare
n }n ∼λ f square ,
with f square (θ1 , θ2 ) = h1D (θ1 )f1D (θ2 ) + f1D (θ1 )h1D (θ2 ), the functions f1D and
h1D being, respectively, the symbols associated to the stiffness and mass matrix
sequences in one dimension, that is, f1D (θ) = 2−2 cos(θ), h1D (θ) = 32 + 13 cos(θ).
Expanding the coefficients, we have
8 2 2 2 2
f square (θ1 , θ2 ) = − cos(θ1 )− cos(θ2 )− cos (θ1 + θ2 )− cos(θ1 −θ2 ). (32)
3 3 3 3 3
Let us now consider {Tn (f square )}n , the sequence of bi-level Toeplitz matrices
generated by f square , which is a multilevel GLT sequence by definition (see [21]).
We adopt the bi-index notation, that is, we consider the elements of Tn (f square )
to be indexed with a bi-index (i1 , i2 ) with i1 = 1, . . . , n1 and i2 = 1, . . . , n2 .
Let us focus on Ωi . In the Q1 case there is a one-to-one correspondence
between the mesh grid points which are in the interior of the subdomain Ωin i
and the Ni degrees of freedom, where Ni is the dimension of Vi,h . See [22,
Section 3] for an explanation.
Following the analysis performed in [5], first, we  set to zero all the rows
(i1 , i2 ) and columns (i1 , i2 ) in Tn (f square ) such that n1i+1
1
, n2i+1
2
̸∈ Ωin
i , which
 
i1 i2
in the Q1 case means that n1 +1 , n2 +1 is not a mesh grid point tied to a
degree of freedom associated with a trial function in Vi,h . We call the matrix
resulting from this operation Ti . Then, we define the restriction maps ΠΩi and
ΠTΩi such that ΠΩi Ti ΠTΩi is the matrix obtained from Ti deleting all rows (i1 , i2 )
 
and columns (i1 , i2 ) such that n1i+1 1
, n2i+1
2
̸∈ Ωin
i . With a similar procedure,
we construct also the matrices ΠΩe Te ΠTΩe . Then, we can apply [5, Lemma 4.4]
and conclude that

{τ ΠΩe Te ΠTΩe }n ∼λ (f e , De ),
{τ ΠΩi Ti ΠTΩi }n ∼λ (f i , Di ),

with f i (x, θ) = τ f square (θ), (x, θ) ∈ Di = Ωi ×[0, π]2 , and f e (x, θ) = τ f square (θ),
(x, θ) ∈ De = Ωe × [0, π]2 . Apart from a low-rank correction due to boundary
conditions, the matrices ΠΩi Ti ΠTΩi coincide with Ain i and the same holds for the
external domain. Since low-rank corrections are zero-distributed and hence in
2 The label Q denotes the space of elements of order p with square support while P is
p p
used for elements in a triangulated mesh. For a broader understanding and the analysis with
the P1 elements, we refer to [5, Section 7].

14
the Hermitian case have no impact on spectral distributions (see Exercise 5.3 in
book [20]), we can safely write

{τ Ain e e
e }n ∼λ (f , D ), (33)
{τ Ain
i }n
i
∼λ (f , D ). i
(34)

Finally, by exploiting again the useful non-uniqueness of the spectral symbol as

in Remark 4.3, we also have

{τ Ain
e }n ∼λ (τ f
square
, [0, π]2 ), (35)
{τ Ain
i }n ∼λ (τ f square 2
, [0, π] ), (36)

since f square over [0, π]2 is a rearrangement of both f e over De and f i over Di .
Given relationships (33) and (34), Theorem 4.3 applies since the other as-
sumptions are verified when using any reasonable discretization. For instance,
the parameter r in Theorem 4.3 is the ratio
µd (Ωi )
, d = 1, 2, 3.
µd (Ωi ) + µd (Ωe )

Furthermore, if µd (Ωi ) = µd (Ωe ), r = 1/2 and Corollary 4.1 applies. We ob-

serve that in our specific setting, due to (35) and (36), also the assumptions
of Corollary 4.2 hold true, so that the assumption on the ratio r is no longer
relevant.

4.4 Monolithic solution strategy

Since, by construction, An is symmetric positive definite (cf. Section 3), we
employ the conjugate gradient method (CG) to solve system (16) iteratively. We
first describe the dependency on τ and subsequently discuss preconditioning.

Conjugate gradient convergence

To better understand the convergence properties of the CG method in the EMI
context, in particular robustness w.r.t. τ , we recall the following result, proved
in [4].

Lemma 4.1 Assume An to be a symmetric positive definite matrix of size n.

Moreover, suppose there exists a, b > 0, both independent of the matrix size n,
and an integer q < n such that

b < λi (An ), for i = 1, ..., q

a ≤ λj (An ) ≤ b, for j = q + 1, ..., n,

then k iterations of CG are sufficient in order to obtain a relative error of ϵ of

the solution, for any initial guess, i.e.

∥ek ∥An /∥e0 ∥An ≤ ϵ,

15
where ek is the error vector at the kth iteration and
k = q + ⌈log(2ϵ−1 )/ log(α−1 )⌉,
√ √
with α = √b−√a .
b+ a

This lemma can be useful if q ≪ n, i.e. if there are q outliers in the spectrum
of An , all larger than b. Interestingly, the convergence of the CG method does
not depend on the magnitude of the outliers, but only on q. The latter is an
explanation of why the convergence rate of the CG applied to a linear system
with matrix τ −1 An does not significantly depend on τ , as we will show in the
numerical section. Indeed, the matrix τ −1 An can be written as
τ −1 An = Bn + R̃n
with
Ain,Γ
   
0 e 0 0 0 0 0 0
AΓ,in
e AΓe 0 0  0 M Γ
0 Γ
Tei
R̃n = τ −1 
 e
Bn =  , ,
 0 0 0 Ain,Γ
i
 0 0 0 0 
Γ T
0 0 AiΓ,in AΓi 0 (Tei ) MiΓ
0
(37)
where the first term does not depend on τ and the second term has rank at
most 2NΓ . Thanks to the Cauchy interlacing theorem (see [11]), if we choose
[a, b] = [λmin (Bn ), λmax (Bn )] in Lemma 4.1, then q ≤ 2NΓ . In practice, for τ
sufficiently small, we have q = 2NΓ , as we show in the numerical section (cf.
Figure 3). Nevertheless, due to the proximity of the lower bound a to zero, the
convergence of the CG method is far from satisfactory, and it becomes essential
to design an appropriate preconditioner for An .

Preconditioning
We propose a theoretical preconditioner, which is proven to be effective through
the spectral analysis carried out in the previous sections. The process of imple-
menting a practical preconditioning strategy will be discussed at the end of the
current section.
Denoting by IeΓ and IiΓ the identity matrices of the same size as AΓe and AΓi
respectively, a first observation is that the block diagonal matrix
 in 
τ Ae 0 0 0
 0 IeΓ 0 0
Pn =  0 in

0 τ Ai 0
0 0 0 IiΓ
is such that the Hermitian matrix
τ Ain,Γ
 
0 e 0 0
τ AΓ,in
e τ A Γ
e + M Γ Γ
e − Ie 0 Γ
Tei 
A n − Pn = 
 0 in,Γ

0 0 τ Ai 
0 Γ T
(Tei ) τ AΓ,in
i
Γ Γ Γ
τ Ai + Mi − Ii

16
has rank less than or equal to 4NΓ = o(n), n = Ni +Ne , with a√similar argument
of Theorem 4.3 proof. In fact the typical behavior is NΓ ∼ n if Ω is a two-
dimensional domain. In the general setting of a d-dimensional domain and with
d−1
any reasonable discretization scheme we have NΓ ∼ n d , d ≥ 1.
The function f square in formula (32) is nonnegative, implying that the ma-
trices Tn (f square ) are Hermitian positive definite according to [21, Theorem 3.1].
The restriction operators ΠΩi and ΠTΩi have full rank and so ΠΩi Tn (f square )ΠTΩi
is Hermitian positive definite. Moreover, in our case the boundary conditions
do not alter the positive definiteness and the same holds for Ain e . Hence, the
−1/2 −1/2
matrix Pn A n Pn is well-defined and we can write

Pn−1/2 An Pn−1/2 = In − Pn−1/2 (Pn − An )Pn−1/2

−1/2 −1/2
with Pn (Pn − An )Pn having rank less than or equal to 4NΓ = o(n),
−1/2 −1/2
n = Ni + Ne , which means that the sequence {Pn (Pn − An )Pn }n is zero-
distributed. Since all the involved matrices are Hermitian, Exercise 5.3 in book
−1/2 −1/2
[20] tells us that the matrix sequence {Pn A n Pn }n is distributed as 1 in
the eigenvalue sense. Therefore Pn could serve as an effective preconditioner for
An , also because, as already pointed out, the number of possible outliers cannot
d−1
grow more than O(n d ), d ≥ 1, d dimensionality of the physical domain.
We now face the task of finding an efficient method for approximating the
matrix-vector products with the inverses of Ain in
e and Ai . Given the multilevel
structure of these matrices and their role as discretization of the Laplacian op-
erator, multigrid techniques are a natural choice. Defining k = ⌊n/4⌋, when
dealing with a bi-level Toeplitz matrix of size n generated by f square , the con-
ventional bisection and linear interpolation operators are an effective choice for
restriction and prolongation operators. These can be paired with the Jacobi
smoother, which requires suitably chosen relaxation parameters. See [30] for
the theoretical background. A similar multigrid strategy can be applied to the
“subdomain” matrices Ain in
e and Ai thanks to [31].
In practice, we will apply one multigrid iteration as preconditioner for the
whole matrix An , avoiding the construction of Pn . It has been proven [31] that
if two linear systems have spectrally equivalent coefficient matrices if a multigrid
method is effective for a linear system, then it is effective also for the other. The
matrix sequences {An }n and {Pn }n share the same spectral distribution, but
the matrices may not be spectrally equivalent due to outliers, which are at most
4NΓ by the Cauchy interlacing theorem. However, outliers are controlled by
the CG method, as shown by Lemma 4.1. For a recent reference on spectrally
equivalent matrix sequences and their use to design fast iterative solvers see [3].

5 Numerical results
5.1 Problem and discretization settings
We consider the following two scenarios:

17
Figure 2: Left: astrocyte geometry and corresponding mesh. A cube with 0.1
sides is shown as a reference. Right: solution corresponding to Table 4.

(i) An idealized geometry for a single cell, i.e. a two-dimensional square

domain Ω = [0, 1]2 with Ωi = [0.25, 0.75]2 , discretized with a uniform
mesh, cf. Figure 1. The domain Ω is discretized considering a uniform
grid with N × N elements. This tessellation results in n = Ni + Ne =
(N p + 1)2 + 2N p degrees of freedom3 , with p the finite element order.
(ii) A three-dimensional geometry, representing one astrocyte [1], cf. Figure 2
discretized with 32365 grid nodes leading to n = 212548 degrees of freedom
for p = 2 (as used in the numerical experiments).

In both cases, we set σi = σe = 1, pure Dirichlet boundary conditions, i.e.

∂ΩD = ∂Ω in (6), and the following right-hand side source in (8):

f (x, y) = sin(2πx) sin(2πy).

5.2 Implementation and solver parameters

We use FEniCS [2, 26] for parallel finite element assembly and multiphenics4 to
handle multiple meshes with common interfaces and the corresponding mixed-
dimensional weak forms. FEniCS wraps PETSc parallel solution strategies. For
multilevel preconditioning, we use hypre algebraic multigrid (AMG) [19], with
default options, cf. Appendix A. For comparative studies, we use MUMPS as a
direct solver and incomplete LU (ILU) preconditioning with zero fill-in. Parallel
ILU relies on a block Jacobi preconditioned, approximating each block inverse
3 The term 2N p = N is present since the membrane degrees of freedom are repeated both
Γ
for Ωi and Ωe .
4 https://multiphenics.github.io/index.html

18
Figure 3: Left: spectrum of τ −1 An varying τ for (N, p) = (16, 1). The
[a, b] = [λmin (Bn ), λmax (Bn )] interval is shown, with reference to Lemma 4.1
and equation (37). Crucially, the majority of eigenvalues are included in [a, b],
which is independent of τ . Right: condition number of An and CG iterations to
convergence. In yellow the theoretical bound is shown, which depends on a, b
and the relative error ϵ.

with ILU. For iterative strategies, we use a tolerance of 10−6 for the relative
unpreconditioned residual, as stopping criteria.
The implementation is verified considering the same benchmark problem as
in [33, Section 3.2], obtaining the expected convergence rates.

5.3 Experiments: eigenvalues and CG convergence

We consider the behaviour of unpreconditioned CG for τ → 0 and scenario (i).
For simplicity, in this section, we use a unit right-hand side, i.e. f = 1/∥1∥2 in
(16). In Figure 3 an example of the spectrum of An is shown; as presented in
Section 4.4, varying τ influences only 2NΓ eigenvalues. Using Lemma 4.1, with
q = 2NΓ , results in a strict upper bound for the number of CG iterations, which
is observed in practice, despite the increasing condition number κ(An ).
In this idealized setting, imposing τ = 1, we numerically confirm the results
of sections 4.2-4.3, showing that the spectral distribution

{An }n ∼λ (g, D),

with the symbol domain

D = [0, 1] × [0, π]2 × [0, π]2

and
g(x, θi , θe ) = f square (θi )ψ[0, 14 ] (x) + f square (θe )ψ( 41 ,1] (x),
reasonably approximates the eigenvalues of An , given a uniform sampling even
for very few grid points in each of the 5 dimensions of D. In order to make

19
 4 4

3 g 3 f square
6j (An ) 6j (An )

2 2

1 1

0 0

0 5000 10000 15000 0 5000 10000 15000

Figure 4: Left: comparison between the eigenvalues λj (An ) and the samples
of g over a uniform grid on D. Right: comparison between the eigenvalues
λj (An ) and the samples of f square over a uniform grid on [0, π]2 . In both cases
(N, p) = (128, 1) and τ = 1.

visualization possible, we evaluate the function g on its domain and then we

arrange the evaluations in ascending order, which corresponds to considering the
monotone rearrangement of g. Note that in this setting the number r defined
in Theorem 4.3 is equal to 1/4. In the left panel in Figure 4 we observe that
the spectrum of An is qualitatively described by the samplings of the function g
over a uniform grid on D. Moreover, in our particular setting, also the spectral
distribution
{An }n ∼λ (f square , [0, π]2 )
holds, and we numerically confirm this result in the right panel of Figure 4,
where the eigenvalues of An are described by the samplings of the function
f square over a uniform grid on [0, π]2 .

5.4 Experiments: preconditioning

In Table 1 we report CG iterations to convergence for scenario (i), varying the
discretization parameters N, p and τ . Similarly, in Table 2, we report iterations
to convergence using the AMG preconditioned CG (PCG), showing robustness
w.r.t. all discretization parameters. We remark that, for PCG, time-to-solution
depends linearly on n, as expected for multilevel methods. In Table 3 we report
runtimes and iterations to convergence for various preconditioning strategies for
scenario (i); assembly and direct solution timings are also reported for practi-
cality. In terms of efficiency, the AMG and ILU preconditioners are preferable.
In Table 4, we show strong scaling data corresponding to scenario (ii), i.e. the
astrocyte geometry. In this case, the ILU preconditioner is preferable both in
terms of runtime and parallel efficiency.
Summing up the numerical results, we can observe how the AMG precondi-
tioned CG exhibit a near-to-optimal robustness, converging in 5 to 7 iterations

20
 (N, p) (32, 1) (64, 1) (128, 1) (256, 1) (512, 1)
n 1153 4353 16897 66561 264193
τ =1 49 95 186 365 718
τ = 10−1 45 90 176 346 686
τ = 10−2 41 79 151 296 583
τ = 10−3 71 122 185 263 487

(N, p) (16, 2) (32, 2) (64, 2) (128, 2) (256, 2)

τ =1 56 110 215 422 830
τ = 10−1 50 104 204 400 793
τ = 10−2 47 92 175 343 674
τ = 10−3 79 129 192 298 563

Table 1: Number of CG iterations to convergence, scenario (i).

for all the cases considered. For time-to-solution, an ILU preconditioning might
be preferable for three-dimensional complex geometries, while AMG remains the
most convenient choice for structured discrete problems. Let us remark that the
ILU preconditioner can be particularly suitable for multiple time steps since the
ILU decomposition has to be computed only once.

6 Concluding remarks
We described numerical approximation schemes for the EMI equations and stud-
ied the structure and spectral features of the coefficient matrices obtained from
a finite element discretization in the so-called single-dimensional formulation.
The obtained spectral information has been employed for designing appropriate
(preconditioned) Krylov solvers. Numerical experiments have been presented
and critically discussed; depending on the context (dimensionality, geometry)
the CG method, preconditioned with AMG or ILU results in an efficient, scal-
able, and robust solution strategy. The spectrum analysis made it possible to
understand the convergence properties, somehow counterintuitive, of CG in this
context.
Due to the intrinsic generality of the tools in Sections 4.1-4.2, the given
machinery can be applied to a variety of discretization schemes including iso-
geometric analysis, finite differences, and finite volumes, in the spirit of the
sections of books [21, 7] dedicated to applications.

Acknowledgments
Pietro Benedusi and Marie Rognes acknowledge support from the Research
Council of Norway via FRIPRO grant #324239 (EMIx) and from the na-
tional infrastructure for computational science in Norway, Sigma2, via grant
#NN8049K. Paola Ferrari and Stefano Serra-Capizzano are partially supported

21
 (N, p) (32, 1) (64, 1) (128, 1) (256, 1) (512, 1)
n 1153 4353 16897 66561 264193
τ =1 5 5 5 6 5
τ = 10−1 5 5 6 6 6
τ = 10−2 5 5 5 6 7
τ = 10−3 5 5 5 6 6

(N, p) (16, 2) (32, 2) (64, 2) (128, 2) (256, 2)

τ =1 5 5 5 6 5
τ = 10−1 5 6 5 6 6
τ = 10−2 4 5 5 5 6
τ = 10−3 6 5 5 6 6

Table 2: Number of PCG iterations to convergence, scenario (i), using a single

AMG iteration as preconditioner.

Runtime [its.]
Assembly 2.6
Direct 1.6
CG 1.4 [583]
PCG(Jacobi) 1.5 [557]
PCG(SOR) 1.2 [199]
PCG(ILU) 0.7 [120]
PCG(AMG) 0.3 [5]

Table 3: Runtimes (s) and iterations to convergence (in square brackets) for
various solution strategies for τ = 0.01 and scenario (i) with (N, p) = (512, 1) for
a single core. For completeness, finite element assembly timing is also reported.

Cores 1 2 4 8
Assembly 8.4 5.5 3.4 2.3
Direct 16 12.1 10.6 9.9
CG 6.3 [725] 3.5 [724] 2.2 [724] 2.0 [720]
PCG(Jacobi) 3.0 [337] 1.7 [337] 1.0 [337] 1.0 [337]
PCG(SOR) 2.6 [129] 2.0 [174] 1.4 [181] 1.2 [185]
PCG(ILU) 2.2 [98] 1.3 [104] 0.8 [113] 0.6 [119]
PCG(AMG) 5.9 [7] 3.9 [7] 3.0 [7] 2.9 [7]

Table 4: Runtimes (s) and iterations to convergence (in square brackets) for
various parallel solution strategies for τ = 0.01 and scenario (ii), cf. the astro-
cyte geometry in Figure 2, and p = 2. For completeness, finite element assembly
timings are also reported.

22
by the Italian Agency INdAM-GNCS. Furthermore, the work of Stefano Serra-
Capizzano is funded from the European High-Performance Computing Joint
Undertaking (JU) under grant agreement No 955701. The JU receives support
from the European Union’s Horizon 2020 research and innovation programme
and Belgium, France, Germany, and Switzerland. Stefano Serra-Capizzano is
also grateful for the support of the Laboratory of Theory, Economics and Sys-
tems – Department of Computer Science at Athens University of Economics
and Business. Special thanks are extended to Jørgen Dokken, Miroslav Kuchta,
Francesco Ballarin, Lars Magnus Valnes, Halvor Herlyng, Abdellah Marwan,
and Marius Causemann for their precious help.

References
[1] Marwan Abdellah, Juan José Garcı́a Cantero, Nadir Román Guerrero,
Alessandro Foni, Jay S Coggan, Corrado Calı̀, Marco Agus, Eleftherios
Zisis, Daniel Keller, Markus Hadwiger, et al. Ultraliser: a framework for
creating multiscale, high-fidelity and geometrically realistic 3D models for
in silico neuroscience. Briefings in bioinformatics, 24(1):bbac491, 2023.
[2] Martin Alnæs, Jan Blechta, Johan Hake, August Johansson, Benjamin
Kehlet, Anders Logg, Chris Richardson, Johannes Ring, Marie E Rognes,
and Garth N Wells. The FEniCS project version 1.5. Archive of Numerical
Software, 3(100):9–23, 2015.

[3] Owe Axelsson, János Karátson, and Frédéric Magoulès. Robust superlinear
krylov convergence for complex noncoercive compact-equivalent operator
preconditioners. SIAM Journal on Numerical Analysis, 61(2):1057–1079,
2023.

[4] Owe Axelsson and Gunhild Lindskog. On the rate of convergence of the pre-
conditioned conjugate gradient method. Numerische Mathematik, 48:499–
523, 1986.
[5] Giovanni Barbarino. A systematic approach to reduced GLT. BIT Numer-
ical Mathematics, 62(3):681–743, 2022.

[6] Giovanni Barbarino, Davide Bianchi, and Carlo Garoni. Constructive ap-
proach to the monotone rearrangement of functions. Expositiones Mathe-
maticae, 40(1):155–175, 2022.
[7] Giovanni Barbarino, Carlo Garoni, and Stefano Serra-Capizzano. Block
generalized locally Toeplitz sequences: theory and applications in the mul-
tidimensional case. Electronic Transactions on Numerical Analysis, 53:113–
216, 2020.
[8] Giovanni Barbarino, Carlo Garoni, and Stefano Serra-Capizzano. Block
generalized locally Toeplitz sequences: theory and applications in the uni-

23
 dimensional case. Electronic Transactions on Numerical Analysis, 53:28–
112, 2020.
[9] Pietro Benedusi, Paola Ferrari, Carlo Garoni, Rolf Krause, and Stefano
Serra-Capizzano. Fast parallel solver for the space-time IgA-DG discretiza-
tion of the diffusion equation. Journal of Scientific Computing, 89(1):20,
2021.
[10] Pietro Benedusi, Carlo Garoni, Rolf Krause, Xiaozhou Li, and Stefano
Serra-Capizzano. Space-time FE-DG discretization of the anisotropic dif-
fusion equation in any dimension: the spectral symbol. SIAM Journal on
Matrix Analysis and Applications, 39(3):1383–1420, 2018.

[11] Rajendra Bhatia. Matrix analysis, volume 169 of Graduate Texts in Math-
ematics. Springer-Verlag, New York, 1997.
[12] Joel A Black, Jeffery D Kocsis, and Stephen G Waxman. Ion channel
organization of the myelinated fiber. Trends in neurosciences, 13(2):48–54,
1990.

[13] Alessio Paolo Buccino, Miroslav Kuchta, Jakob Schreiner, and Kent-André
Mardal. Improving Neural Simulations with the EMI Model, pages 87–98.
Springer, 2021.
[14] Giacomo Rosilho de Souza, Rolf Krause, and Simone Pezzuto. Boundary
integral formulation of the cell-by-cell model of cardiac electrophysiology.
arXiv preprint arXiv:2302.05281, 2023.
[15] Ali Dorostkar, Maya Neytcheva, and Stefano Serra-Capizzano. Spectral
analysis of coupled PDEs and of their Schur complements via generalized
locally Toeplitz sequences in 2D. Computer Methods in Applied Mechanics
and Engineering, 309:74–105, 2016.
[16] S-E Ekström and Stefano Serra-Capizzano. Eigenvalues and eigenvectors
of banded Toeplitz matrices and the related symbols. Numerical Linear
Algebra with Applications, 25(5):e2137, 2018.
[17] Ada J Ellingsrud, Cécile Daversin-Catty, and Marie E Rognes. A cell-based
model for ionic electrodiffusion in excitable tissue, pages 14–27. Springer
International Publishing, 2021.
[18] Ada J Ellingsrud, Andreas Solbrå, Gaute T Einevoll, Geir Halnes, and
Marie E Rognes. Finite element simulation of ionic electrodiffusion in cel-
lular geometries. Frontiers in Neuroinformatics, 14:11, 2020.

[19] Robert D Falgout and Ulrike Meier Yang. hypre: A library of high per-
formance preconditioners. In Computational Science—ICCS 2002: Inter-
national Conference Amsterdam, The Netherlands, April 21–24, 2002 Pro-
ceedings, Part III, pages 632–641. Springer, 2002.

24
[20] Carlo Garoni, Stefano Serra-Capizzano, et al. Generalized locally Toeplitz
sequences: theory and applications, volume 1. Springer, 2017.
[21] Carlo Garoni, Stefano Serra-Capizzano, et al. Generalized locally Toeplitz
sequences: theory and applications, volume 2. Springer, 2018.

[22] Carlo Garoni, Stefano Serra-Capizzano, and Debora Sesana. Spectral anal-
ysis and spectral symbol of d-variate Qp Lagrangian FEM stiffness matri-
ces. SIAM Journal on Matrix Analysis and Applications, 36(3):1100–1128,
2015.
[23] Ngoc Mai Monica Huynh, Fatemeh Chegini, Luca Franco Pavarino, Martin
Weiser, and Simone Scacchi. Convergence analysis of BDDCprecondition-
ers for hybrid DG discretizations of the cardiac cell-by-cell model. arXiv
preprint arXiv:2212.12295, 2022.
[24] Karoline Horgmo Jæger, Andrew G Edwards, Wayne R Giles, and Aslak
Tveito. From millimeters to micrometers; re-introducing myocytes in mod-
els of cardiac electrophysiology. Frontiers in Physiology, 12:763584, 2021.
[25] Karoline Horgmo Jæger, Kristian Gregorius Hustad, Xing Cai, and Aslak
Tveito. Efficient numerical solution of the EMI model representing the
extracellular space (E), cell membrane (M) and intracellular space (I) of a
collection of cardiac cells. Frontiers in Physics, 8:579461, 2021.

[26] Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution
of differential equations by the finite element method: The FEniCS book,
volume 84. Springer Science & Business Media, 2012.
[27] Mariarosa Mazza, Ahmed Ratnani, and Stefano Serra-Capizzano. Spectral
analysis and spectral symbol for the 2D curl-curl (stabilized) operator with
applications to the related iterative solutions. Mathematics of Computation,
88(317):1155–1188, 2019.
[28] Mariarosa Mazza, Matteo Semplice, Stefano Serra-Capizzano, and Elena
Travaglia. A matrix-theoretic spectral analysis of incompressible Navier–
Stokes staggered DG approximations and a related spectrally based pre-
conditioning approach. Numerische Mathematik, 149(4):933–971, 2021.
[29] Yoichiro Mori and Charles Peskin. A numerical method for cellular electro-
physiology based on the electrodiffusion equations with internal boundary
conditions at membranes. Communications in Applied Mathematics and
Computational Science, 4(1):85–134, 2009.

[30] Stefano Serra-Capizzano and Cristina Tablino-Possio. Multigrid methods

for multilevel circulant matrices. SIAM Journal on Scientific Computing,
26(1):55–85, 2004.

25
[31] Stefano Serra-Capizzano and Cristina Tablino-Possio. Two-grid methods
for Hermitian positive definite linear systems connected with an order re-
lation. Calcolo, 51(2):261–285, 2014.
[32] Paolo Tilli. A note on the spectral distribution of Toeplitz matrices. Linear
and Multilinear Algebra, 45(2-3):147–159, 1998.
[33] Aslak Tveito, Karoline H Jæger, Miroslav Kuchta, Kent-Andre Mardal,
and Marie E Rognes. A cell-based framework for numerical modeling of
electrical conduction in cardiac tissue. Frontiers in Physics, page 48, 2017.
[34] Aslak Tveito, Kent-André Mardal, and Marie E Rognes. Modeling excitable
tissue: the EMI framework. Springer, 2021.
[35] Evgenij E Tyrtyshnikov. A unifying approach to some old and new the-
orems on distribution and clustering. Linear algebra and its applications,
232:1–43, 1996.

A Appendix: multigrid solver parameters

We provide the default AMG parameters in PETSc obtained using the ksp view
command.

KSP Object: 1 MPI process

type: cg
maximum iterations=10000, initial guess is zero
tolerances: relative=1e-06, absolute=1e-50, divergence=10000.
left preconditioning
using UNPRECONDITIONED norm type for convergence test
PC Object: 1 MPI process
type: hypre
HYPRE BoomerAMG preconditioning
Cycle type V
Maximum number of levels 25
Maximum number of iterations PER hypre call 1
Convergence tolerance PER hypre call 0.
Threshold for strong coupling 0.25
Interpolation truncation factor 0.
Interpolation: max elements per row 0
Number of levels of aggressive coarsening 0
Number of paths for aggressive coarsening 1
Maximum row sums 0.9
Sweeps down 1
Sweeps up 1
Sweeps on coarse 1
Relax down symmetric-SOR/Jacobi
Relax up symmetric-SOR/Jacobi

26
 Relax on coarse Gaussian-elimination
Relax weight (all) 1.
Outer relax weight (all) 1.
Using CF-relaxation
Not using more complex smoothers.
Measure type local
Coarsen type Falgout
Interpolation type classical
SpGEMM type cusparse
linear system matrix = precond matrix

27