Abstract
In this work, we are interested in solving large linear systems stemming from the extra–membrane–intra model, which is employed for simulating excitable tissues at a cellular scale. After setting the related systems of partial differential equations equipped with proper boundary conditions, we provide its finite element discretization 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 spectral information is used for designing appropriate preconditioned Krylov solvers. Through numerical experiments, we show that the presented solution strategy is robust w.r.t. problem and discretization parameters, efficient and scalable.
Similar content being viewed by others
Avoid common mistakes on your manuscript.
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 cardiology and neuroscience, where spreading excitation by electrical signalling plays a crucial role [13, 16, 17, 24, 25, 29, 32, 34]. We refer to [35] 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.
The latter information is crucial to understand the algebraic properties of such discrete objects and to design iterative solution strategies, enabling the efficient solution of large scale models, as in [9, 10]. 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 Sect. 2 the continuous problem is introduced together with possible generalizations. Section 3 considers a basic Galerkin strategy for approximating the examined problem. The spectral analysis is given in Sect. 4 in terms of distribution results and degenerating eigenspaces. Specifically, Sect. 4.1 lays out the foundational theories and concepts necessary for understanding the distribution of the spectrum of the entire system (presented in a general form in Sect. 4.2) and the specific stiffness matrices and matrix sequences (discussed in Sect. 4.3). These findings are then employed in Sect. 4.4 for proposing a preconditioning strategy. In Sect. 5 we report and critically discuss the numerical experiments and Sect. 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 \(\varGamma \) 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 \(\varOmega =\varOmega _i\cup \varOmega _e\cup \varGamma \subset \mathbb {R}^d\), typically with \(d\in \{2,3\}\), \(\partial \varOmega =\partial \varOmega _e/\varGamma \), and \(\partial \varOmega _i=\varGamma \), we consider the following stationary problem for the intra- and extra-cellular potentials \(u_i,u_e\), and for the membrane current \(I_m\):
with \(\sigma _e,\sigma _i,\tau \in \mathbb {R}_+\), \(\textbf{n}_i\) (resp. \(\textbf{n}_e\)) is the outer normal on \(\partial \varOmega _i\) (resp. \(\partial \varOmega _e\)) and \(f\in L^2(\varGamma )\) is known. We can close the EMI problem with homogeneous boundary conditions:
with \(\partial \varOmega =\partial \varOmega _D\cup \partial \varOmega _N\). In the case of a pure Neumann problem, with \(\partial \varOmega =\partial \varOmega _N\), uniqueness must be enforced through an additional constraint (e.g. a condition on the integral of \(u_e\)). We are imposing homogeneous boundary conditions 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 \(\varGamma \) through a Robin-type condition (3)–(5), depending on the current \(I_m\). The EMI problem can be considered a mixed-dimensional problem since the unknowns of interest are defined on sets of different dimensionality, d for \(u_i\) and \(u_e\) and \(d-1\) for \(I_m\).
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\in (0,T]\) given a final time \(T>0\):
where the positive constant \(C_m\) is a capacitance and the ionic current \(I_{\textrm{ion}}\) is a reaction term. An implicit (resp. explicit) integration of \(I_m\) (resp. \(I_{\textrm{ion}}\)), with time step \(\varDelta t>0\), results in Eq. (5) with
and \(\tau = C_m^{-1}\varDelta t\).
The model in Eqs. (1)–(5) can describe a single cell or multiple disjoint cells, i.e. with \(\varOmega _i=\bigcup _{j=1}^{N_\text {cell}}\varOmega _i^j\), in an extracellular media, cf. Fig. 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 [35] 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 \(\varOmega _i,\varOmega _e\), and possibly \(\varGamma \).
We use the so-called single-dimensional formulation and the corresponding discrete operators. In this setting, the weak form depends only on bulk quantities \(u_i\) and \(u_e\) since the current term \(I_m\) is replaced by:
according to equations (4)–(5). Let us remark that “single” refers to the previous substitution, eliminating the variable defined in \(\varGamma \); the overall EMI problem is still in multiple dimensions, in the sense that \(d\ge 1\).
After substituting the expression for \(I_m\) in (3), assuming the solution \(u_r\), for \(r\in \{i,e\}\), to be sufficiently regular over \(\varOmega _r\), we multiply the PDEs in (1)–(2) by test functions \(v_r(\textbf{x})\in V_r(\varOmega _r)\), with \(V_r\) a sufficiently regular Hilbert space with elements satisfying the boundary conditions in (6)–(7); in practice \(H^1(\varOmega _i)\) and \(H^1(\varOmega _e)\), would be a standard choice. After integrating over \(\varOmega _r\) and applying integration by parts, using the normal flux definition (3), the weak EMI problem reads: find \(u_i\in V_i(\varOmega _i)\) and \(u_e\in V_e(\varOmega _e)\) such that
for all test functions \(v_e\in V_e(\varOmega _e)\) and \(v_i\in V_i(\varOmega _i)\). We refer to [35, Section 6.2.1] for boundedness and coercivity results for this formulation.
For each subdomain \(\varOmega _r\), with \(r\in \{i,e\}\), we construct a conforming tassellation \(\mathcal {T}_r\). We then introduce a yet unspecified discretization via finite element basis functions (e.g. Lagrangian elements of order \(p\in \mathbb {N}\) on a regular grid) for \(V_e\) and \(V_i\):
with \(N_e,N_i\in \mathbb {N}\) denoting the number of degrees of freedom in the corresponding subdomains. We can further decompose \(N_r=N_{r,\text {in}}+N_{r,\varGamma }\), i.e. further differentiating between internal and membrane degrees of freedom, with \(N_{r,\varGamma }\) basis functions \(\phi ^r_j\) having support intersecting \(\varGamma \), cf. Fig. 1. In the numerical experiments, we consider matching \(\mathcal {T}_r\) on the interface \(\varGamma \) and the same p for \(V_{e,h}\) and \(V_{i,h}\); nevertheless Sects. 3 and 4 are developed in a general setting, so that the theory is ready also for potential extensions.
Left: example of EMI geometry with \(\varOmega _i=\varOmega _i^1\cup \varOmega _i^2\). Right: example of the discretization setting for \(d=2\) on a regular grid. We notice that interior and exterior nodes overlap in \(\varGamma \). For \(p = 1\), in this case, we have \(N_i = 9, N_e = 24\) with \(N_{i,\text {in}}= 1, N_{e,\text {in}} = 16\), and \(N_\varGamma =8\). Corresponding labels are reported
From (9)–(10) we define the following discrete operators: intra- and extra- Laplacians
membrane mass matrices:
and the coupling matrix
Finally, we write the linear system of size \(n = N_e+N_i\) corresponding to (9)–(10) as
with \(\textbf{u}_r,\textbf{f}_r\in \mathbb {R}^{N_r}\) the unknowns and the right hand side corresponding to \(\varOmega _r\) and \(r\in \{i,e\}\). We remark that the operator \(A_n\) is symmetric and positive definite.Footnote 1 Making the degrees of freedom in \(\varGamma \) denoted explicitly (with reference to Fig. 1), i.e. \(\textbf{u}_{i,\varGamma },\textbf{u}_{e,\varGamma }\in \mathbb {R}^{N_\varGamma }\), as well as the interior ones, i.e. \(\textbf{u}_{i,\textrm{in}}\in \mathbb {R}^{N_i - N_\varGamma }\) and \(\textbf{u}_{e,\textrm{in}}\in \mathbb {R}^{N_e - N_\varGamma }\), we can rewrite more extensively system (16) as
with
For forthcoming use, we define the bulk matrix
where all membrane terms are zeroed, and the difference matrix
which contains all the terms corresponding to the membrane \(\varGamma \).
4 Spectral Analysis
In this section, we study the spectral distribution of the matrix sequence \(\{A_{n}\}_n\) under various assumptions for determining the global behaviour of the eigenvalues of \(A_n\) 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 [7, 8, 19, 20].
First, we give the formal definition of Toeplitz structures, eigenvalue (spectral) 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 provide 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,
As indicated in (21), left, the matrix T can be infinite-dimensional. A finite-dimensional Toeplitz matrix of dimension n is denoted as \(T_n\). We also consider sequences of Toeplitz matrices as a function of their dimension, denoted as \(\{T_n\}_n,\, n=1,2,\ldots ,\infty \).
In general, the entries \(a_k\), \(k\in {\mathbb Z}\), can be matrices (blocks) themselves, defining \(T_n\) as a block Toeplitz matrix. Thus, in the block case, \(a_k, k\in \{-n+1, \ldots , n-1\}\) are blocks of size \(s_1\times s_2\), the subindex of \(T_n\) is the number of blocks in the Toeplitz matrix while \(N_1\times N_2\) is the size of the matrix with \(N_1 = n\,s_1\), \(N_2 = n\,s_2\). To be specific, we use also the notation \(X_{N_1,N_2}=T_n\). A special case of block-Toeplitz matrices is the class of two- and multilevel block Toeplitz matrices, 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 \(Q^d=[-\pi ,\pi ]^d, d\ge 1\), with d-dimensional Lebesgue measure \(\mu _d (Q^d)=(2\pi )^d\). Denote by \(f_k\) the Fourier coefficients of f,
where \(\theta =(\theta _1,\ldots ,\theta _d)\), \((k,{ \theta })=\sum _{j=1}^{d}k_j \theta _j\), \(n=(n_1,\ldots ,n_d)\), and \(N(n)={n_1}\times \cdots \times n_d\). By following the multi-index notation in [36, Section 6], with each f we can associate a sequence of Toeplitz matrices \(\{T_n\}_n\), where
\(\textbf{e}=[1,1,\ldots ,1]\in {\mathbb N}^d\). For \(d=1\)
or for \(d = 2\), i.e. the two-level case, and for example \(n=(2,3)\), we have
The function f is referred to as the generating function (or the symbol of) \(T_n\). Using a more compact notation, we say that the function f is the generating function of the whole sequence \(\{T_n\}_n\) and we write \(T_n = T_n(f)\).
If f is d-variate, \(\mathbb {C}^{s_1\times s_2}\) matrix-valued, and integrable over \(Q^d\), \(d,s_1,s_2\ge 1\), i.e. \(f\in L^1(Q^d,s_1\times s_2)\), then we can define the Fourier coefficients of f in the same way (now \(f_k\) is a matrix of size \(s_1\times s_2\)) and consequently \(T_n = \{ f_{k-\ell }\}_{k,\ell =\textbf{e}^T}^{n} \in \mathbb {C}^{s_1N(n)\times s_2N(n)}\), then \(T_n\) is a d-level block Toeplitz matrix according to Definition 4.1. If \(s_1=s_2=s\) then we write \(f\in L^1(Q^d, s)\).
As in the scalar case, the function f is referred to as the generating function of \(T_n\). We say that the function f is the generating function of the whole sequence \(\{T_n\}_n\), and we use the notation \(T_n = T_n(f)\).
Definition 4.3
Let \(f:D\rightarrow \mathbb {C}^{s\times s}\) be a measurable matrix-valued function with eigenvalues \(\lambda _i(f)\) and singular values \(\sigma _i(f)\), \(i=1,\ldots ,s\). Assume that \(D\subset \mathbb {R}^d\) is Lebesgue measurable with positive and finite Lebesgue measure \(\mu _d(D)\). Assume that \(\{A_n\}_n\) is a sequence of matrices such that \(\textrm{dim}(A_n)=d_n\rightarrow \infty \), as \(n\rightarrow \infty \) and with eigenvalues \(\lambda _j(A_n)\) and singular values \(\sigma _j(A_n)\), \(j=1,\ldots ,d_n\).
-
We say that \(\{A_n\}_{n}\) is distributed as f over D in the sense of the eigenvalues, and we write \(\{A_n\}_{n}\sim _\lambda (f,D),\) if
$$\begin{aligned} \lim _{n\rightarrow \infty }\frac{1}{d_n}\sum _{j=1}^{d_n}F(\lambda _j(A_n))= \frac{1}{\mu _d(D)} \int _D \frac{1}{s}\sum _{i=1}^sF(\lambda _i(f(t))) \, \textrm{d}t, \end{aligned}$$(22)for every continuous function F with compact support. In this case, we say that f is the spectral symbol of \(\{A_{n}\}_{n}\).
-
We say that \(\{A_n\}_{n}\) is distributed as f over D in the sense of the singular values, and we write \(\{A_n\}_{n}\sim _\sigma (f,D)\), if
$$\begin{aligned} \lim _{n\rightarrow \infty }\frac{1}{d_n}\sum _{j=1}^{d_n}F(\sigma _j(A_n))= \frac{1}{\mu _d(D)} \int _D \frac{1}{s}\sum _{i=1}^sF(\sigma _i(f(t))) \, \textrm{d}t, \end{aligned}$$(23)for every continuous function F with compact support. In this case, we say that f is the singular value symbol of \(\{A_{n}\}_{n}\).
-
The notion \(\{A_n\}_{n}\sim _\sigma (f,D)\) applies also in the rectangular case where f is \(\mathbb {C}^{s_1\times s_2}\) matrix-valued. In such a case the parameter s in formula (23) has to be replaced by the minimum between \(s_1\) and \(s_2\): furthermore \(A_n\in \mathbb {C}^{d_n^{(1)}\times d_n^{(2)}}\) with \(d_n\) in formula (23) being the minimum between \(d_n^{(1)}\) and \(d_n^{(2)}\). 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 \(\{A_n\}_n\sim _{\lambda ,\sigma }(f,D)\) with \(\{A_n\}_n\sim _{\lambda ,\sigma } f\).
Remark 4.1
If f is smooth enough, an informal interpretation of the limit relation (22) (resp. (23)) is that when n is sufficiently large, the eigenvalues (resp. singular values) of \(A_{n}\) can be subdivided into s different subsets of the same cardinality. Then \(d_n/s\) eigenvalues (resp. singular values) of \(A_{n}\) can be approximated by a sampling of \(\lambda _1(f)\) (resp. \(\sigma _1(f)\)) on a uniform equispaced grid of the domain D, and so on until the last \(d_n/s\) eigenvalues (resp. singular values), which can be approximated by an equispaced sampling of \(\lambda _s(f)\) (resp. \(\sigma _s(f)\)) in the domain D.
Remark 4.2
We say that \(\{A_n\}_n\) is zero-distributed in the sense of the eigenvalues if \(\{A_n\}_n\sim _\lambda 0\). Of course, if the eigenvalues of \(A_n\) tend all to zero for \(n\rightarrow \infty \), then this is sufficient to claim that \(\{A_n\}_n\sim _\lambda 0\).
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 [7, 20] and references therein).
Theorem 4.1
[33] Let \(f\in L^1(Q^d,s_1\times s_2)\), then \(\{T_{n}(f)\}_{{n}}\sim _\sigma (f,Q^d).\) If \(s_1=s_2=s\) and if f is a Hermitian matrix-valued function, then \(\{T_{n}(f)\}_{{n}}\sim _\lambda (f,Q^d)\).
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 \(\{A_n\}_n\) be a sequence of matrices, with \(A_n\) Hermitian of size \(d_n\), and let \(\{P_n\}_n\) be a sequence such that \(P_n\in \mathbb C^{d_n\times \delta _n}\), \(P_n^*P_n=I_{\delta _n}\), \(\delta _n\le d_n\) and \(\delta _n/d_n\rightarrow 1\) as \(n\rightarrow \infty \). Then \(\{A_n\}_n\sim _{\lambda }f\) if and only if \(\{P_n^*A_nP_n\}_n\sim _{\lambda }f\).
With the notations of the result above, the matrix sequence \(\{P_n^*A_nP_n\}_n\) is called a compression of \(\{A_n\}_n\) and the single matrix \(P_n^*A_nP_n\) is called a compression of \(A_n\).
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 [14] and the discussion below Theorem 3 at p. 8 in [15].
Remark 4.3
[Multiple Toeplitz generating functions and multiple spectral symbols] Let n be a multiple of k and consider the Toeplitz matrix \(X_n = T_n(f)\). We can view \(X_n\) as a block-Toeplitz with blocks of size \(k \times k\) such that \(X_n = T_n(f) = T_{\frac{n}{k}}(f^{[k]})\). In our specific context, we have
with \({\textbf{e}}_j\), \(j=1,\ldots ,k\), being the canonical basis of \(\mathbb {C}^k\).
As an example, let n be even and consider the function \(f(\theta )=2-2\cos (\theta )\). According to Definition 4.2, \(X_n=T_n(f)\) but we can also view the matrix as \(X_n = T_{\frac{n}{2}}(f^{[2]})\), namely
where
Thus, \(f^{[2]} = \left( A_0+A_1 e^{i\theta } +A_1^T e^{-i\theta }\right) \).
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 \(\{T_{n}(f)\}_{{n}}\sim _\lambda (f,Q^d)\) and \(\{T_{n}(f)\}_{{n}}\sim _\lambda (f^{[k]},Q^d)\) 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 [6, 15]). A simple case is given by standard Toeplitz sequences \(\{T_n(f)\}_n\), with f real-valued and even that is \(f(\theta )=f(-\theta )\) almost everywhere, \(\theta \in Q\). In that case
due to the even character of f, and hence it is also true that \(\{T_n(f)\}_n\sim _{\lambda } (f,Q_+)\), \(Q_+=(0,\pi )\). A general analysis of these concepts via rearrangement theory can be found in [6] and references therein.
Remark 4.4
We remark that the presented tools are general and can be applied to matrix sequences stemming from a variety of discretization schemes such as isogeometric analysis or finite volumes, in the spirit of the sections of books [7, 20] dedicated to applications and of the exposition paper [22].
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 \(N_i,N_e\) and \(N_\varGamma \) denote the number of intra-, extra- and membrane degrees of freedom. We recall that the system matrix \(A_n\) and its decomposition as \({\tilde{A}}_n+R_n\) are defined in Eqs. (18)–(20).
Theorem 4.3
Assume that
Assume that
\(D^e\subset \mathbb {R}^{k_e}\), \(D^i\subset \mathbb {R}^{k_i}\), given \(f^e\) and \(f^i\) extra- and intra- spectral symbols. It follows that
where \(g(x,t_i,t_e)=f^i(t_i) \psi _{[0,r]}(x) + f^e(t_e) \psi _{(r,1]}(x)\), \(x\in [0,1]\), \(t_i\in D^i\), \(t_e\in D^e\), and
\(\psi _Z\) denoting the characteristic function of the set Z.
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 generality we assume that both \(f^e\), \(f^i\) take values into \(\mathbb {C}^{s\times s}\) (the case where \(f^e\), \(f^i\) have different matrix sizes is treated in Remark 4.6, for the sake of notational simplicity). Because of Eqs. (25) and (26), setting \(N_e'=N_e-N_\varGamma \), \(N_i'=N_i-N_\varGamma \), and taking into account Definition 4.3, we have
for any continuous function F with bounded support. Now we build the \(2\times 2\) block matrix
and we consider again a generic continuous function F with bounded support. By the block diagonal structure of \(X[N_e,N_i,N_\varGamma ]\), we infer
As a consequence, taking the limit as \(N_i,N_e\rightarrow \infty \), using the assumption \(N_\varGamma = o(\min \{N_i,N_e\})\),
we obtain that the limit of \(\varDelta (N_e,N_i,N_\varGamma ,F)\) exists and it is equal to
With regard to Definition 4.3, the quantity (29) does not look like the right-hand side of (22), since we see two different terms. The difficulty can be overcome by enlarging the space with a fictitious domain [0, 1] and interpreting the sum in (29) as the global integral involving a step function. In reality, we rewrite (29) as
with \(D=[0,1]\times D^i\times D^e\), \(q=k_e+k_i+1\), \(\mu _{q}(D)=\mu _{k_e}(D^e)\mu _{k_i}(D^i)\). Hence the proof of the relation
is complete. Now we observe that \(X[N_e,N_i,N_\varGamma ]\) is a compression of \({\tilde{A}}_n\) (cf. definition after Theorem 4.2): therefore by exploiting again the hypothesis \(N_\varGamma = o(\min \{N_i,N_e\})\) which implies \(N_\varGamma = o(n)\), by Theorem 4.2, we deduce
Finally, in [19] the following is proven: if \(\{V_n\}_n\sim _{\lambda }f\) (not necessarily Toeplitz or GLT), \(\{A_n=V_n+Y_n\}_n\) with \(\{Y_n\}_n\) zero-distributed and \(A_n,Y_n\) Hermitian for every dimension, then \(\{A_n\}_n\sim _{\lambda }f\). This setting is exactly our setting since by direct inspection the rank of
is bounded by \(4N_\varGamma \). Hence the number of nonzero eigenvalues of the latter matrix cannot exceed \(4N_\varGamma =o(n)\), \(n=N_i+N_e\), and therefore the related matrix sequence is zero-distributed in the eigenvalue sense i.e. \(\{R_n\}_n\sim _\lambda 0\). Since \(A_n={\tilde{A}}_n + R_n\), the application of the statement in [19, Exercise 5.3] implies directly \(\{A_n\}_n\sim _\lambda (g,[0,1]\times D^i\times D^e)\) and the proof is concluded. \(\square \)
Remark 4.5
An informal implication of Theorem 4.3 is that, under the listed assumptions, the spectrum of \(A_n\) is essentially driven by the bulk blocks, while membrane terms are negligible for large n, i.e. \(\{R_n\}_n \sim _\lambda 0.\)
Remark 4.6
The case where \(f^e\) and \(f^i\) have different matrix sizes would complicate 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 [6, 15] 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
Assume that
\(D\subset \mathbb {R}^{k}\). It follows that
and
Proof
Notice that the writing
and
are equivalent for \(c=1/2\). \(\square \)
Corollary 4.2
Assume that
Assume that
\(D\subset \mathbb {R}^{k}\). It follows that
and
Proof
The proof of the relation \(\{{\tilde{A}}_n\}_n,\ \ \{A_n\}_n\sim _\lambda (f,D)\) follows directly from the limit displayed in (29), after replacing \(f^e,f^i\) with f and necessarily \(D^e,D^i\) with D. The rest is obtained verbatim as in Theorem 4.3. \(\square \)
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 preconditioning approach. According to [21, Section 5], given a bi-dimensional square domain, the sequence of corresponding stiffness matrices \(\{A_n^\mathrm{{square}}\}_n\) obtained from the linear \(\mathbb {Q}_1\) Lagrangian finite element methodFootnote 2 (FEM) has the following spectral distribution:
with \(f^\mathrm{{square}}(\theta _1,\theta _2) = h_\mathrm{{1D}}(\theta _1) f_\mathrm{{1D}}(\theta _2) + f_\mathrm{{1D}}(\theta _1) h_\mathrm{{1D}}(\theta _2)\), the functions \(f_\mathrm{{1D}}\) and \(h_\mathrm{{1D}}\) being, respectively, the symbols associated to the stiffness and mass matrix sequences in one dimension, that is, \(f_\mathrm{{1D}}(\theta ) = 2-2\cos (\theta )\), \(h_\mathrm{{1D}}(\theta ) = \frac{2}{3}+\frac{1}{3}\cos (\theta )\). Expanding the coefficients, we have
Let us now consider \(\{T_n(f^\mathrm{{square}})\}_n\), the sequence of bi-level Toeplitz matrices generated by \(f^\mathrm{{square}}\), which is a multilevel GLT sequence by definition (see [20]). We adopt the bi-index notation, that is, we consider the elements of \(T_n(f^\mathrm{{square}})\) to be indexed with a bi-index \((i_1,i_2)\) with \(i_1 = 1,\dots ,n_1\) and \(i_2 = 1,\dots ,n_2\).
Let us focus on \(\varOmega _i\). In the \(\mathbb {Q}_1\) case there is a one-to-one correspondence between the mesh grid points which are in the interior of the subdomain \(\varOmega _i^\textrm{in}\) and the \(N_i\) degrees of freedom, where \(N_i\) is the dimension of \(V_{i,h}\). See [21, Section 3] for an explanation.
Following the analysis performed in [5], first, we set to zero all the rows \((i_1,i_2)\) and columns \((i_1,i_2)\) in \(T_n(f^\mathrm{{square}})\) such that \(\left( \frac{i_1}{n_1+1},\frac{i_2}{n_2+1}\right) \not \in \varOmega _i^\textrm{in}\), which in the \(\mathbb {Q}_1\) case means that \(\left( \frac{i_1}{n_1+1},\frac{i_2}{n_2+1}\right) \) is not a mesh grid point tied to a degree of freedom associated with a trial function in \(V_{i,h}\). We call the matrix resulting from this operation \(T_i\). Then, we define the restriction maps \(\varPi _{\varOmega _i}\) and \(\varPi _{\varOmega _i}^T\) such that \(\varPi _{\varOmega _i} T_i \varPi _{\varOmega _i}^T\) is the matrix obtained from \(T_i\) deleting all rows \((i_1,i_2)\) and columns \((i_1,i_2)\) such that \(\left( \frac{i_1}{n_1+1},\frac{i_2}{n_2+1}\right) \not \in \varOmega _i^\textrm{in}\). With a similar procedure, we construct also the matrices \(\varPi _{\varOmega _e} T_e \varPi _{\varOmega _e}^T\). Then, we can apply [5, Lemma 4.4] and conclude that
with \(f^i(\textbf{x},\varvec{\theta }) = \tau f^\mathrm{{square}}(\varvec{\theta })\), \((\textbf{x},\varvec{\theta }) \in D^i = \varOmega _i \times [0,\pi ]^2\), and \(f^e(\textbf{x},\varvec{\theta }) = \tau f^\mathrm{{square}}(\varvec{\theta })\), \((\textbf{x},\varvec{\theta }) \in D^e = \varOmega _e \times [0,\pi ]^2\). Apart from a low-rank correction due to boundary conditions, the matrices \(\varPi _{\varOmega _i} T_i \varPi _{\varOmega _i}^T\) coincide with \(A_{i}^{\textrm{in}}\) and the same holds for the external domain. Since low-rank corrections are zero-distributed and hence in the Hermitian case have no impact on spectral distributions (see Exercise 5.3 in book [19]), we can safely write
Finally, by exploiting again the useful non-uniqueness of the spectral symbol as in Remark 4.3, we also have
since \(f^\mathrm{{square}}\) over \([0,\pi ]^2\) is a rearrangement of both \(f^e\) over \(D^e\) and \(f^i\) over \(D^i\).
Given relationships (35) and (36), Theorem 4.3 applies since the other assumptions are verified when using any reasonable discretization. For instance, the parameter r in Theorem 4.3 is the ratio
Furthermore, if \(\mu _d(\varOmega _i)=\mu _d(\varOmega _e)\), \(r=1/2\) and Corollary 4.1 applies. We observe that in our specific setting, due to (37) and (38), 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, \(A_n\) is symmetric positive definite (cf. Sect. 3), we employ the conjugate gradient method (CG) to solve system (16) iteratively. We first describe the dependency on \(\tau \) and subsequently discuss preconditioning.
4.4.1 Conjugate Gradient Convergence
To better understand the convergence properties of the CG method in the EMI context, in particular robustness w.r.t. \(\tau \), we recall the following result, proved in [4].
Lemma 4.1
Assume \(A_n\) 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
then k iterations of CG are sufficient in order to obtain a relative error of \(\epsilon \) of the solution, for any initial guess, i.e.
where \(\textbf{e}^k\) is the error vector at the kth iteration and
with \(\alpha =\frac{\sqrt{b}-\sqrt{a}}{\sqrt{b}+\sqrt{a}}\).
This lemma can be useful if \(q\ll n\), i.e. if there are q outliers in the spectrum of \(A_n\), 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 \(\tau ^{-1}A_n\) does not significantly depend on \(\tau \), as we will show in the numerical section. Indeed, the matrix \(\tau ^{-1}A_n\) can be written as
with
where the first term does not depend on \(\tau \) and the second term has rank at most \(2N_\varGamma \). Thanks to the Cauchy interlacing theorem (see [11]), if we choose \([a,b]=[\lambda _{\min }(B_n),\lambda _{\max }(B_n)]\) in Lemma 4.1, then \(q \le 2N_\varGamma \). In practice, for \(\tau \) sufficiently small, we have \(q = 2N_\varGamma \), as we show in the numerical section (cf. Fig. 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 \(A_n\).
4.4.2 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 implementing a practical preconditioning strategy is discussed at the end of the current section.
Denoting by \(I_{e}^{\varGamma }\) and \(I_{i}^{\varGamma }\) the identity matrices of the same size as \(A_{e}^{\varGamma }\) and \(A_{i}^{\varGamma }\) respectively, a first observation is that the block diagonal matrix
is such that the Hermitian matrix
has rank less than or equal to \(4N_\varGamma =o(n)\), \(n=N_i+N_e\), using an argument similar to the one in the proof of Theorem 4.3. In fact the typical behavior is \(N_\varGamma \sim \sqrt{n}\) if \(\varOmega \) is a two-dimensional domain. In the general setting of a d-dimensional domain and with any reasonable discretization scheme we have \(N_\varGamma \sim n^{d-1\over d}\), \(d\ge 1\).
The function \(f^\mathrm{{square}}\) in formula (34) is nonnegative, implying that the matrices \(T_n(f^\mathrm{{square}})\) are Hermitian positive definite according to [20, Theorem 3.1]. The restriction operators \(\varPi _{\varOmega _i}\) and \(\varPi _{\varOmega _i}^T\) have full rank and so \(\varPi _{\varOmega _i}T_n(f^\mathrm{{square}})\varPi _{\varOmega _i}^T\) is Hermitian positive definite. Moreover, in our case the boundary conditions do not alter the positive definiteness and the same holds for \(A_{e}^{\textrm{in}}\). Hence, the matrix \(P_n^{-1/2}A_nP_n^{-1/2}\) is well-defined and we can write
with \(P_n^{-1/2}(P_n-A_n)P_n^{-1/2}\) having rank less than or equal to \(4N_\varGamma =o(n)\), \(n=N_i+N_e\), which means that the sequence \(\{P_n^{-1/2}(P_n-A_n)P_n^{-1/2}\}_n\) is zero-distributed. Since all the involved matrices are Hermitian, Exercise 5.3 in book [19] tells us that the matrix sequence \(\{P_n^{-1/2}A_nP_n^{-1/2}\}_n\) is distributed as 1 in the eigenvalue sense. Therefore \({P_n}\) could serve as an effective preconditioner for \({A_n}\), also because, as already pointed out, the number of possible outliers cannot grow more than \(O(n^{d-1\over d})\), \(d\ge 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 \(A_{e}^{\textrm{in}}\) and \(A_{i}^{\textrm{in}}\). Given the multilevel structure of these matrices and their role as discretization of the Laplacian operator, multigrid techniques are a natural choice. Defining \(k=\lfloor n/4 \rfloor \), when dealing with a bi-level Toeplitz matrix of size n generated by \(f^\mathrm{{square}}\), the conventional 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 \(A_{e}^{\textrm{in}}\) and \(A_{i}^{\textrm{in}}\) thanks to [31].
In practice, we will apply one multigrid iteration as preconditioner for the whole matrix \(A_n\), avoiding the construction of \(P_n\). It has been proven [31] that if two linear systems have coefficient matrices that are spectrally equivalent, the effectiveness of a multigrid method for one system also extends to the other. The matrix sequences \(\{A_n\}_n\) and \(\{P_n\}_n\) share the same spectral distribution, but the matrices may not be spectrally equivalent due to outliers, which are at most \(4N_\varGamma \) 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:
-
(i)
An idealized geometry for a single cell, i.e. a two-dimensional square domain \(\varOmega =[0,1]^2\) with \(\varOmega _i=[0.25,0.75]^2\), discretized with a uniform mesh, cf. Fig. 1. The domain \(\varOmega \) is discretized considering a uniform grid with \(N\times N\) elements. This tessellation results in \(n=N_i+N_e=(Np+1)^2+2Np\) degrees of freedom,Footnote 3 with p the finite element order.
-
(ii)
A three-dimensional geometry, representing one astrocyte [1], cf. Fig. 2 discretized with 32,365 grid nodes leading to \(n=212548\) degrees of freedom for \(p = 2\) (as used in the numerical experiments).
In both cases, we set \(\sigma _i=\sigma _e=1\), pure Dirichlet boundary conditions, i.e. \(\partial \varOmega _D=\partial \varOmega \) in (6), and the following right-hand side source in (8):
Left: astrocyte geometry and corresponding mesh. A cube with 0.1 sides is shown as a reference. Right: solution corresponding to Table 4
5.2 Implementation and Parameters for the Solution Methods
We use FEniCS [2, 26] for parallel finite element assembly and multiphenicsFootnote 4 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 boomerAMG algebraic multigrid (AMG) [18], with default options, cf. Appendix A. In the 3D case, i.e. scenario (ii), for AMG we use a threshold of 0.9 for the strong coupling.Footnote 5 since since the default value (0.25) is optimized for 2D applications. 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 with ILU. In all cases, preconditioners are directly applied to the system matrix \(A_n\). 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 [34, Section 3.2], obtaining the expected convergence rates.
5.3 Experiments: Eigenvalues and CG Convergence
We consider the behaviour of unpreconditioned CG for \(\tau \rightarrow 0\) and scenario (i). For simplicity, in this section, we use a unit right-hand side, i.e. \(\textbf{f}=\textbf{1}/\Vert \textbf{1}\Vert _2\) in (16). In Fig. 3 an example of the spectrum of \(A_n\) is shown; as presented in Sect. 4.4, varying \(\tau \) influences only \(2N_\varGamma \) eigenvalues. Using Lemma 4.1, with \(q = 2N_\varGamma \), results in a strict upper bound for the number of CG iterations, which is observed in practice, despite the increasing condition number \(\kappa (A_n)\).
Left: spectrum of \(\tau ^{-1}A_n\) varying \(\tau \) for \((N,p)=(16,1)\). The \([a,b]=[\lambda _{\min }(B_n), \lambda _{\max }(B_n)]\) interval is shown, with reference to Lemma 4.1 and Eq. (39). Crucially, the majority of eigenvalues are included in [a, b], which is independent of \(\tau \). Right: condition number of \(A_n\) and CG iterations to convergence. In yellow the theoretical bound is shown, which depends on a, b and the relative error \(\epsilon \)
In this idealized setting, imposing \(\tau = 1\), we numerically confirm the results of Sects. 4.2–4.3, showing that the spectral distribution
with the symbol domain
and
reasonably approximates the eigenvalues of \(A_n\), given a uniform sampling even for very few grid points in each of the 5 dimensions of D. In order to make 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 Fig. 4 we observe that the spectrum of \(A_{n}\) 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
holds, and we numerically confirm this result in the right panel of Fig. 4, where the eigenvalues of \(A_{n}\) are described by the samplings of the function \(f^\mathrm{{square}}\) over a uniform grid on \([0,\pi ]^2\).
Left: comparison between the eigenvalues \(\lambda _j(A_{n})\) and the samples of g over a uniform grid on D. Right: comparison between the eigenvalues \(\lambda _j(A_{n})\) and the samples of \(f^\mathrm{{square}}\) over a uniform grid on \([0,\pi ]^2\). In both cases \((N,p)=(128,1)\) and \(\tau = 1\)
5.4 Experiments: Preconditioning
In Fig. 5, we report the spectra of \(A_n\) and \(P_n^{-1}A_n\) to highlight the clustering of eigenvalues around one in the preconditioned case. Such behavior is more precisely quantified in Table 1, where, for \(\epsilon >0\), we report the number of outliers given by \(c_n=c_n(\epsilon )\), the latter being the number of eigenvalues of the preconditioned matrix not belonging to \([1-\epsilon ,1+\epsilon ]\). As expected, we observe \(c_n/n \rightarrow 0\) as the total size n is increased and \(c_n\sim \sqrt{n}\), showing the effectiveness of the theoretical preconditioner \(P_n\).
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 practicality. In terms of efficiency, the AMG preconditioner is preferable. In Table 4, we show strong scaling data corresponding to scenario (ii), i.e. the astrocyte geometry. In this case, AMG and ILU preconditioners are both competitive in terms of runtime and parallel efficiency.
Summing up the numerical results, we can observe how the AMG-preconditioned CG exhibits a near-to-optimal robustness, converging in 5 to 7 iterations in 2D and in 12 iterations in 3D. Regarding time-to-solution, AMG is also the most convenient choice, with ILU being a close competitor. Regarding scalability, a parallel efficiency of approximately 80% is reached for both approaches, before saturation.
6 Concluding Remarks
We described numerical approximation schemes for the EMI equations and studied 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; the CG method, preconditioned with AMG results in an efficient, scalable, and robust solution strategy. The spectrum analysis made it possible to understand the convergence properties, somehow counterintuitive, of CG in this context.
Data Availability
Data sets generated during the current study are available from the corresponding author on reasonable request.
Notes
Dirichlet boundary conditions can be imposed to the right-hand side to enforce symmetry.
The label \(\mathbb {Q}_p\) denotes the space of elements of order p with square support while \(\mathbb {P}_p\) is used for elements in a triangulated mesh. For a broader understanding and the analysis with the \(\mathbb {P}_1\) elements, we refer to [5, Section 7].
The term \(2Np=N_\varGamma \) is present since the membrane degrees of freedom are repeated both for \(\varOmega _i\) and \(\varOmega _e\).
With the command -pc_hypre_boomeramg_strong_threshold.
References
Abdellah, M., Cantero, J.J.G., Guerrero, N.R., Foni, A., Coggan, J.S., Calì, C., Agus, M., Zisis, E., Keller, D., Hadwiger, M., et al.: Ultraliser: a framework for creating multiscale, high-fidelity and geometrically realistic 3D models for in silico neuroscience. Brief. Bioinform. 24(1), bbac491 (2023)
Alnæs, M., Blechta, J., Hake, J., Johansson, A., Kehlet, B., Logg, A., Richardson, C., Ring, J., Rognes, M.E., Wells, G.N.: The FEniCS project version 1.5. Arch. Numer. Softw. 3(100), 9–23 (2015)
Axelsson, O., Karátson, J., Magoulès, F.: Robust superlinear Krylov convergence for complex noncoercive compact-equivalent operator preconditioners. SIAM J. Numer. Anal. 61(2), 1057–1079 (2023)
Axelsson, O., Lindskog, G.: On the rate of convergence of the preconditioned conjugate gradient method. Numer. Math. 48, 499–523 (1986)
Barbarino, G.: A systematic approach to reduced GLT. BIT Numer. Math. 62(3), 681–743 (2022)
Barbarino, G., Bianchi, D., Garoni, C.: Constructive approach to the monotone rearrangement of functions. Expos. Math. 40(1), 155–175 (2022)
Barbarino, G., Garoni, C., Serra-Capizzano, S.: Block generalized locally Toeplitz sequences: theory and applications in the multidimensional case. Electron. Trans. Numer. Anal. 53, 113–216 (2020)
Barbarino, G., Garoni, C., Serra-Capizzano, S.: Block generalized locally Toeplitz sequences: theory and applications in the unidimensional case. Electron. Trans. Numer. Anal. 53, 28–112 (2020)
Benedusi, P., Ferrari, P., Garoni, C., Krause, R., Serra-Capizzano, S.: Fast parallel solver for the space–time IgA-DG discretization of the diffusion equation. J. Sci. Comput. 89(1), 20 (2021)
Benedusi, P., Garoni, C., Krause, R., Li, X., Serra-Capizzano, S.: Space–time FE-DG discretization of the anisotropic diffusion equation in any dimension: the spectral symbol. SIAM J. Matrix Anal. Appl. 39(3), 1383–1420 (2018)
Bhatia, R.: Matrix Analysis. Graduate Texts in Mathematics, vol. 169. Springer, New York (1997). https://doi.org/10.1007/978-1-4612-0653-8
Black, J.A., Kocsis, J.D., Waxman, S.G.: Ion channel organization of the myelinated fiber. Trends Neurosci. 13(2), 48–54 (1990)
Buccino, A.P., Kuchta, M., Schreiner, J., Mardal, K.A.: Improving Neural Simulations with the EMI Model, pp. 87–98. Springer, Berlin (2021)
Dorostkar, A., Neytcheva, M., Serra-Capizzano, S.: Spectral analysis of coupled PDEs and of their Schur complements via generalized locally Toeplitz sequences in 2D. Comput. Methods Appl. Mech. Eng. 309, 74–105 (2016)
Ekström, S.E., Serra-Capizzano, S.: Eigenvalues and eigenvectors of banded Toeplitz matrices and the related symbols. Numer. Linear Algebra Appl. 25(5), e2137 (2018)
Ellingsrud, A.J., Daversin-Catty, C., Rognes, M.E.: A Cell-Based Model for Ionic Electrodiffusion in Excitable Tissue. Springer, Berlin (2021)
Ellingsrud, A.J., Solbrå, A., Einevoll, G.T., Halnes, G., Rognes, M.E.: Finite element simulation of ionic electrodiffusion in cellular geometries. Front. Neuroinform. 14, 11 (2020)
Falgout, R.D., Yang, U.M.: Hypre: a library of high performance preconditioners. In: Computational Science-ICCS 2002: International Conference Amsterdam, The Netherlands, April 21–24, 2002 Proceedings, Part III, pp. 632–641. Springer (2002)
Garoni, C., Serra-Capizzano, S.: Generalized Locally Toeplitz Sequences: Theory and Applications, vol. 1. Springer, Berlin (2017)
Garoni, C., Serra-Capizzano, S.: Generalized Locally Toeplitz Sequences: Theory and Applications, vol. 2. Springer, Berlin (2018)
Garoni, C., Serra-Capizzano, S., Sesana, D.: Spectral analysis and spectral symbol of \(d\)-variate \(\mathbb{Q} _{\textbf{p} }\) Lagrangian FEM stiffness matrices. SIAM J. Matrix Anal. Appl. 36(3), 1100–1128 (2015)
Garoni, C., Speleers, H., Ekström, S.E., Reali, A., Serra-Capizzano, S., Hughes, T.J.R.: Symbol-based analysis of finite element and isogeometric B-spline discretizations of eigenvalue problems: exposition and review. Arch. Comput. Methods Eng. 26(5), 1639–1690 (2019)
Huynh, N.M.M., Chegini, F., Pavarino, L.F., Weiser, M., Scacchi, S.: Convergence analysis of BDDCpreconditioners for hybrid DG discretizations of the cardiac cell-by-cell model. arXiv preprint arXiv:2212.12295 (2022)
Jæger, K.H., Edwards, A.G., Giles, W.R., Tveito, A.: From millimeters to micrometers; re-introducing myocytes in models of cardiac electrophysiology. Front. Physiol. 12, 763584 (2021)
Jæger, K.H., Hustad, K.G., Cai, X., Tveito, A.: 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. Front. Phys. 8, 579461 (2021)
Logg, A., Mardal, K.A., Wells, G.: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, vol. 84. Springer, Berlin (2012)
Mazza, M., Ratnani, A., Serra-Capizzano, S.: Spectral analysis and spectral symbol for the 2D curl–curl (stabilized) operator with applications to the related iterative solutions. Math. Comput. 88(317), 1155–1188 (2019)
Mazza, M., Semplice, M., Serra-Capizzano, S., Travaglia, E.: A matrix-theoretic spectral analysis of incompressible Navier–Stokes staggered DG approximations and a related spectrally based preconditioning approach. Numer. Math. 149(4), 933–971 (2021)
Mori, Y., Peskin, C.: A numerical method for cellular electrophysiology based on the electrodiffusion equations with internal boundary conditions at membranes. Commun. Appl. Math. Comput. Sci. 4(1), 85–134 (2009)
Serra-Capizzano, S., Tablino-Possio, C.: Multigrid methods for multilevel circulant matrices. SIAM J. Sci. Comput. 26(1), 55–85 (2004)
Serra-Capizzano, S., Tablino-Possio, C.: Two-grid methods for Hermitian positive definite linear systems connected with an order relation. Calcolo 51(2), 261–285 (2014). https://doi.org/10.1007/s10092-013-0081-9
de Souza, G.R., Krause, R., Pezzuto, S.: Boundary integral formulation of the cell-by-cell model of cardiac electrophysiology. arXiv preprint arXiv:2302.05281 (2023)
Tilli, P.: A note on the spectral distribution of Toeplitz matrices. Linear Multilinear Algebra 45(2–3), 147–159 (1998)
Tveito, A., Jæger, K.H., Kuchta, M., Mardal, K.A., Rognes, M.E.: A cell-based framework for numerical modeling of electrical conduction in cardiac tissue. Front. Phys. 5, 48 (2017)
Tveito, A., Mardal, K.A., Rognes, M.E.: Modeling Excitable Tissue: The EMI Framework. Springer, Berlin (2021)
Tyrtyshnikov, E.E.: A unifying approach to some old and new theorems on distribution and clustering. Linear Algebra Appl. 232, 1–43 (1996)
Acknowledgements
Pietro Benedusi and Marie Rognes acknowledge support from the Research Council of Norway via FRIPRO Grant #324239 (EMIx) and from the national infrastructure for computational science in Norway, Sigma2, via Grant #NN8049K. Paola Ferrari and Stefano Serra-Capizzano are partially supported 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 Systems - Department of Computer Science at Athens University of Economics and Business. Special thanks are extended to the anonymous Referee for the enriching comments and to Jørgen Dokken, Miroslav Kuchta, Francesco Ballarin, Lars Magnus Valnes, Halvor Herlyng, Abdellah Marwan, and Marius Causemann for their precious help.
Funding
Open access funding provided by Università della Svizzera italiana
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors declare that they have no conflict of interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix: Parameters of Multigrid
Appendix: Parameters of Multigrid
We provide the default AMG parameters in PETSc obtained using the ksp_view command. For the 3D case, i.e. scenario (ii), the threshold for strong coupling is set to 0.9.
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
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
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Benedusi, P., Ferrari, P., Rognes, M.E. et al. Modeling Excitable Cells with the EMI Equations: Spectral Analysis and Iterative Solution Strategy. J Sci Comput 98, 58 (2024). https://doi.org/10.1007/s10915-023-02449-2
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10915-023-02449-2