In the EMI model for a single cell, like illustrated in Figure 1, we consider an intracellular domain, Ω i , surrounded by an extracellular domain, Ω e . The boundary between the intracellular and extracellular domains defines the cell membrane, Γ . The EMI model (see, e.g., Ref. 8) describes the electric potential in these domains and is given by the following system of equations: ∇ ⋅ σ e ∇ u e = 0     in   Ω e , (1) ∇ ⋅ σ i ∇ u i = 0     in   Ω i , (2) u e = 0     at   ∂ Ω e D , (3) n e ⋅ σ e ∇ u e = 0     at   ∂ Ω e N , (4) n e ⋅ σ e ∇ u e = − n i ⋅ σ i ∇ u i ≡ I m   at   Γ , (5) u i − u e = v     at   Γ , (6) v t = 1 C m ( I m − I ion )   at   Γ , (7) s t = F     at   Γ . (8) Here, u e is the extracellular potential defined in Ω e , u i is the intracellular potential defined in Ω i , and v is the membrane potential defined at Γ . The potentials are all given in units of mV. The parameters σ e and σ i represent the extracellular and intracellular conductivites, respectively, and are given in units of mS/cm. Furthermore, n e and n i denote the outward pointing unit normal vectors of the extracellular and intracellular spaces, respectively, and C m represents the specific membrane capacitance, given in units of μF/cm². Time is given in units of ms and length is given in units of cm.

The current density I ion represents the sum of the ionic current densities through different types of ion channels, pumps and exchangers on the cell membrane, and I m represents the sum I ion and the capacitive current density C m v t . All current densities are given in units of μA/cm². The model for I ion typically involves a set of additional state variables describing the gating mechanisms of the ion channels and ionic concentrations. These variables are defined on Γ and are denoted by s in the EMI system defined above. The dynamics of s are governed by a set of ordinary differential equations described by F ( v , s ) . In the computations of this paper, we let I ion and F ( v , s ) be given by the Grandi et al. model for human ventricular cardiomyocytes [33].

On the outer boundary of the extracellular space, we apply the homogeneous Dirichlet boundary condition Eq. 3 on a part of the boundary, ∂ Ω e D , and the homogeneous Neumann boundary condition Eq. 4 on the remaining part, ∂ Ω e N . In the computations reported below, we apply Dirichlet boundary conditions on the outer extracellular boundary in the x- and y-directions and Neumann boundary conditions on the extracellular boundary in the z-direction.

In order to illustrate how the linear part of the EMI model Eqs 1–8 can be discretized and solved in a coupled manner, we describe below how to set up a finite difference discretization of the full EMI system. However, because of the possibly nonlinear terms included in F and I ion , we use a temporal operator splitting method to split the EMI system into a linear part and a nonlinear part. The discrete version of this approach is taken from Ref. 8 but repeated here for completeness.

The coupled discrete version of the EMI model may be used to find accurate solutions of the system (see, e.g., Refs 8–10). However, in some cases it would be much more convenient to be able to rewrite the EMI system as a set of more standard problems that can be solved individually using existing well-developed solution methods for such standard problems. Linear systems arising from discretizations of elliptic partial differential equations have been under intense scrutiny for decades and the efforts have resulted in optimal solution strategies; see, e.g., Ref. 30. Our aim is therefore to introduce a spatial operator splitting algorithm for the EMI model that results in sub-problems formulated in terms of standard elliptic problems, thus enabling the use of optimal solvers.

In this section, we describe a method for constructing such a splitting of the EMI model. In this approach, for each time step, we will solve the equations of the intracellular and extracellular spaces separately as standard Laplace equations with Robin, Neumann or Dirichlet boundary conditions on the boundary of the domains. The algorithm described in this section is summarized in Algorithm 1.

An illustration of the computational mesh used in each of the steps of the finite difference discretization of the splitting method is found in Figure 3B. In the discrete version of the splitting approach, we solve Step 1 just like for the coupled version by taking m forward Euler steps of size Δ t * = Δ t / m . In Step 2, we consider the nodes representing the membrane and the purely intracellular points and we have one equation and one unknown, u ¯ i , in all these points. For the purely intracellular points, we discretize Eq. 17 using Eq. 11. For the membrane points, we discretize n i ⋅ σ i ∇ u ¯ i in Eq. 18 using a standard first-order finite difference approximation. In Step 3, we discretize Eqs 20–22 in the same manner as for the coupled system, and use standard first-order finite difference approximations of n i ⋅ σ i ∇ u ¯ i and n e ⋅ σ e ∇ u ¯ e in Eqs 19 and 23, respectively. In Step 4, we use the straightforward discrete version of Eq. 26.

In order to investigate the accuracy of the splitting method introduced above, we set up a simple test case with a single cell surrounded by an extracellular space. The cell is located in the center of the domain, and the size of the entire domain Ω = Ω i ∪ Ω e is 180 μm × 130 μm × 26 μm. The default parameter values used in the simulation are given in Table 1, and the membrane dynamics are modeled by the epicardial version of the Grandi et al. model for human ventricular cardiomyocytes [33]. In addition, 20% of the cell in the lower x-direction is stimulated by a membrane current of −40 μA/cm².

Figures 4 and 5 show a comparison of the solutions of the coupled approach and the splitting approach applied to this simple test problem. The left panel shows the extracellular potential, u e , and the right panel shows the intracellular potential, u i . Furthermore, the upper panel shows the solution of the coupled approach, the middle panel shows the solution of the splitting approach, and the lower panel shows the absolute value of the difference between the two solutions. In Figure 4 we use one iteration ( N it = 1 ) in the splitting approach, and in Figure 5, we use five iterations ( N it = 5 ) . We observe that the coupled and splitting approaches result in very similar solutions for both u e and u i and for both one and five iterations of the splitting approach. However, the difference is reduced when five iterations are used compared to when one iteration is used. For N it = 1 , the maximum difference between the two solutions is about 10 − 6  mV for u e and 10 − 2  mV for u i . For N it = 5 , the maximum difference is reduced to about 10 − 14  mV for u e and 10 − 10  mV for u i .

In Table 2, we report the difference between the solutions of the two approaches in the form of the relative, discrete l 2 -norm of the difference. We consider a number of different values of Δ t for the splitting approach and compare these solutions to the solution of the coupled approach for Δ t = 0.001  ms. In the splitting approach, we take five iterations for each time step ( N it = 5 ) . In the table, we observe that as Δ t is reduced for the splitting approach, the difference between the two solutions decreases with a rate roughly equal to Δ t .

A two-dimensional version of the domain for the EMI model for two connected cells is illustrated in Figure 6. Two cells, Ω i 1 and Ω i 2 , are surrounded by an extracellular space, Ω e , and the boundary between Ω e and each of the cells Ω i 1 and Ω i 2 consists of the cell membranes denoted by Γ 1 and Γ 2 , respectively. In addition, the boundary between the two intracellular domains defines an intercalated disc, denoted by Γ 1,2 . The EMI model describes the electric potential in such a domain and is given by the equations ∇ ⋅ σ e ∇ u e = 0   in   Ω e , (27) ∇ ⋅ σ i ∇ u i k = 0   in   Ω i k , (28) u e = 0   at   ∂ Ω e D , (29) n e ⋅ σ e ∇ u e = 0   at   ∂ Ω e N , (30) n e ⋅ σ e ∇ u e = − n i k ⋅ σ i ∇ u i k ≡ I m k   at   Γ k , (31) u i k − u e = v k   at   Γ k , (32) v t k = 1 C m ( I m k − I ion k )   at   Γ k , (33) s t k = F k   at   Γ k , (34) n i k ˜ ⋅ σ i ∇ u i k ˜ = − n i k ⋅ σ i ∇ u i k ≡ I k , k ˜   at   Γ k , k ˜ , (35) u i k − u i k ˜ = w k   at   Γ k , k ˜ , (36) w t k = 1 C g ( I k , k ˜ − I gap k )   at   Γ k , k ˜ , (37) for each cell k and each neighboring cell k ˜ . Here, u e , u i k , v k are the extracellular, intracellular and membrane potentials in Ω e , Ω i k and at Γ k , respectively. In addition, w k = u i k − u i k ˜ denotes the potential difference at an intercalated disc, Γ k , k ˜ . Note here that at a given intercalated disc, this potential difference is defined differently for the two adjoining cells. For instance, at the intercalated disc Γ 1,2 illustrated in Figure 6, w 1 = u i 1 − u i 2 and w 2 = u i 2 − u i 1 . The current density I gap k (in μA/ cm 2 ) represents the ionic current density through the gap junctions at the intercalated disc and is given by I gap k = 1 R g w k , (38) where R g is the gap junction resistance (in k Ω cm 2 ), and C g in Eq. 37 represents the specific capacitance of the intercalated disc. Furthermore, I k , k ˜ (in μA/ cm 2 ) represents the sum of I gap k and the capacitive current over the gap junction, C g w t k . For definitions and units for the remaining variables and parameters, see Section 2.1.

A two-dimensional version of the computational mesh used in the finite difference discretization of the EMI model for two connected cells is illustrated in Figure 7A. We observe that in this case, there are four different types of nodes; nodes located in the extracellular space, marked by × , nodes located in the intracellular space, marked by ○, nodes located at the membrane, marked by ⊗ , and, finally, nodes located at the intercalated disc, marked by • .

In the coupled discrete version of the EMI model for connected cells, we split the EMI model into a nonlinear part and a linear part using the same temporal operator splitting procedure as for the discrete coupled model for a single cell (see Section 2.2). Moreover, the first nonlinear step, solving the ordinary differential equations describing the membrane dynamics, is solved in exactly the same manner as for the single cell case (see Section 2.2.2).

Furthermore, in the finite difference discretization of the linear part of the EMI model we use the exact same approach to discretize the equations for the intracellular, extracellular and membrane nodes as for the single cell case (see Section 2.2.3). The only exceptions requiring special attention in the case of connected cells are the extracellular nodes located directly below or above the intercalated disc (see Figure 7A). For these nodes, we enforce a no-flow boundary condition and set the extracellular potential equal to the potential in the purely extracellular nodes located directly above or below the corner nodes.

The nodes located at the intercalated disc are treated in exactly the same manner as the membrane nodes except that I gap k is not set to zero in Eq. 37. We discretize n i k ˜ ⋅ σ i ∇ u i k ˜ and n i k ⋅ σ i ∇ u i k in Eq. 35 using standard first-order finite differences, and use a straightforward discrete version of Eq. 36. Finally, Eq. 37 is discretized in time using an implicit finite difference approximation w k , n − w k , n − 1 Δ t = 1 C g ( I k , k ˜ n − 1 R g w k , n ) , where a first-order finite difference approximation of − n i k ⋅ σ i ∇ u i k , n is used to define I k , k ˜ n .

We now want to split the intracellular subdomains in order to obtain classical elliptic equations defined on small domains. In this splitting approach, the linear part of the EMI model is for each time step split into standard Laplace problems with Robin, Neumann or Dirichlet boundary conditions. These problems can be solved separately in each of the domains Ω i k and Ω e .

Steps 1, 3, and 4 of the splitting approach are exactly the same for connected cells as for the single cell case (see Section 2.3), but Step 2 is extended to account for the cell coupling at the intercalated discs, Γ k , k ˜ . In this section we therefore only describe the extended version of Step 2 and refer to Section 2.3 for descriptions of the remaining steps of the procedure. In addition, the full splitting algorithm for coupled cells is summarized in Algorithm 2.

Figure 7B shows an illustration of the computational mesh used in each of the steps of the finite difference discretization of the splitting approach for connected cells. The discrete versions of steps 1, 3, and 4 of the procedure are set up in the same manner as for a single cell (see Section 2.4). For each cell and in each iteration of Step 2, we consider the nodes representing the intercalated disc, the membrane and the purely intracellular points, and we have one equation and one unknown, u ¯ i k , in all these nodes. For the purely intracellular nodes, we discretize Eq. 39 using Eq. 11. For the membrane and intercalated disc nodes, we discretize n i ⋅ σ i ∇ u ¯ i k in Eqs 40 and 41 using a standard first-order finite difference.

In order to investigate the accuracy of the splitting method for the EMI model for connected cells, we extend the simple test case introduced in Section 2.5 to include a collection of 5 × 5 connected cells. We still use the parameters specified in Table 1, unless otherwise specified. In order to allow for two-dimensional cell connections, each cell is extended with additional cubes of smaller width in the x- and y-directions (see Figure 8). The sizes of these additional cubes are given in Table 3. The distance between the intracellular domain and the extracellular boundary is 12 μm in the x- and y-directions and 4 μm in the z-direction. Furthermore, we stimulate the 2 × 2 cells in the lower left corner by a 1 ms-long stimulus current of − 40  μA/ cm 2 , and we consider the solution at t = 2.5  ms.

Note that in the current implementation of the model, the intracellular and extracellular conductivities, σ i and σ e , are assumed to be isotropic, and the anisotropy of cardiac tissue is represented by the anisotropic cell shape and distribution of gap junctions. Anisotropic conductivities, σ i and σ e , could, however, also be straightforwardly incorporated in the model (see Eqs 27–37). Using the parameters of the simple test case setup, the conduction velocity is approximately 26 cm/s in the longitudinal direction (x-direction) and approximately 6 cm/s in the transverse direction (y-direction).

In Figure 9, we compare the solutions of the coupled and splitting approaches applied to this test problem. In the splitting approach, we use five inner and outer iterations ( M i t = N i t = 5 ) . We observe that the coupled and splitting approaches result in very similar solutions for both u e and u i . The maximum difference between the two solutions is about 10 − 8  mV for u e and 10 − 6  mV for u i .

Table 4 compares the solution of the splitting approach for different values of Δ t to the solution of the coupled approach using Δ t = 0.001  ms. Again, we use five inner and outer iterations ( M i t = N i t = 5 ) in the splitting approach. We observe that as the size of Δ t is decreased for the splitting approach, the difference between the two solutions decreases with a rate roughly equal to Δ t . Note that in this case, Δ t = 0.1  ms lead to instabilities for the splitting approach, so Δ t = 0.05  ms is the smallest time step considered.

In Figure 10, we investigate how the absolute value of the difference between the intracellular potentials, computed using the two solution methods, depends on the time step, Δ t , the number of iterations used for the coupling between the intracellular and extracellular domains, N it , and the number of iterations used for the coupling between the interior of the cells, M it . In this case, we use the same Δ t in the coupled and splitting approaches. In the first row of Figure 10, we observe that when only one inner and outer iteration is used in the splitting approach, the maximum difference between the solutions is in the range of 1–10 mV, and that the difference decreases as the time step is decreased. In the next rows, we observe that the difference between the solutions of the two methods decreases as we increase the number of iterations used in the splitting approach. For Δ t = 0.05  ms, we see that the maximum difference is decreased from about 8 mV to about 1.5 mV if N it or M it is increased to five. Moreover, if both N it and M it is set to five, the maximum difference decreases to less than 0.01 mV. For Δ t = 0.01  ms and Δ t = 0.001  ms, we observe that the difference decreases more extensively when the number of outer iterations, N it , is increased to five than when the number of inner iterations, M it , is increased to five. Furthermore, when both are set to five the maximum difference is about 10 − 6 for Δ t = 0.01  ms and 10 − 4 for Δ t = 0.001  ms.

Figure 11 shows a similar comparison of the solutions of the membrane potential, v, at the center of the membrane of the center cell as a function of time. We consider three different time steps Δ t = 0.05  ms, Δ t = 0.01  ms and Δ t = 0.001  ms and some different combinations of the number of inner and outer iterations in the splitting algorithm. The upper panel of the figure shows the full upstroke of the action potential. In order to get a better impression of the differences between the solutions, the lower panel zooms in on the solution at the action potential peak. In the left panel, we observe that when only one iteration is used in the splitting approach ( N i t = M i t = 1 ) , there are clearly visible differences between the solutions of the splitting approach and the coupled approach. However, in the next panels, we observe that when the number of iterations is increased, the difference between the solutions decreases, and for N it = M it = 5 , the solutions of the coupled and splitting approaches are completely indistinguishable in the plots.

In Figure 12, we examine how the maximum difference in the computed intracellular potentials between the coupled and splitting approaches decreases as the number of iterations used in the splitting approach increases. In the left panel, we fix the number of inner iterations, M it , at one or five and vary the number of outer iterations, N it , between one and eight. We observe that the difference between the solutions of the two methods decreases as the number of outer iterations, N it , increases. In addition, the difference is smaller if the number of inner iterations, M it , is five than if it is one. In the right panel, we fix the number of outer iterations, N it , at one or five and vary the number of inner iterations, M it , between one and eight. In this case we observe that increasing M it above two or three when N it is fixed at one or five, respectively, does not seem to increase the accuracy of the splitting approach.

Recall that the intracellular space consists of the individual excitable cells. Due to the splitting algorithm, the intracellular systems associated with the cells can be solved independently. This immediately brings a cell-level parallelism, because the individual intracellular systems can be assigned to the different CPU cores. With a spatial resolution of Δ x = Δ y = Δ z = 2   μ m there are about 5,300 degrees of freedom in the linear system for computing the intracellular potential of each cell. Such a small-size system suits well for a direct solver, see, e.g., Ref. 36, for which the LU factorization is an excellent choice. Moreover, since each intracellular matrix A i that arises from the finite difference discretization remains unchanged throughout an entire EMI simulation, the factorization step is executed only once. The resulting lower and upper triangular matrices, L i and U i , are stored in memory, so that they are used in a forward-backward substitution procedure every time a linear system of the form A i x = b associated with an intracellular subdomain needs to be solved.

It is possible to adopt several threads to parallelize each forward-backward substitution step, but the amount of parallelism at this level is limited. Such a parallelization approach is thus not recommended unless the number of available CPU cores greatly exceeds the number of cells.

For computing the LU factorization and performing forward-back substitutions, we have chosen the serial version of the SuperLU library [37]. We assign each of the linear systems A i x = b to an OpenMP thread so that multiple threads can work in parallel solving different linear systems. SuperLU was configured to use the included BLAS implementation, CBLAS, for the linear algebra routines. It is important that the BLAS implementation used by SuperLU is “thread safe” (i.e., that it supports being called from multiple threads) and does not attempt any parallelization of its own, as this interferes with the cell-level parallelization strategy employed here.

Apart from the aforementioned advantage that the LU factorization of each A i is only computed once, there is another opportunity for performance enhancement. This is true as long as we are interested in studying excitable cells of the same shape and size. Under such an assumption, the number of unique intracellular matrices is limited. More specifically, if the excitable cells form a rectangular 2D grid, then for each of the four sides, a cell can either be connected to a neighboring cell via a gap junction or not, leading to some slight changes in the cell’s geometry and boundary condition. As such, there are at most 2 4 = 16 unique matrices, and in the examples in this paper where the cells form an M × N rectangular grid with M , N ≥ 3 , there are only 9 different matrices (one for each of the four corners, one for each of the four sides, and one for the inner cells). To save memory usage, we can thus group together the linear systems sharing the same matrix A i , and instead solve A i X = B , where X and B are matrices that contain, respectively, all the relevant solution vectors and right-hand sides. Furthermore, we only need to compute and store the LU factorization for each of the 9 unique matrices. In order to ensure good data locality, the CPU cores that are assigned to the same intracellular matrix A i should preferably have direct access to the memory domain where L i and U i reside.

One remark is in order here. Although we discuss in this paper a two dimensional collection of cells, it is important to keep in mind that all the cell geometries are in 3D and all cells are surrounded by a 3D extracellular space. In fact, the EMI model for a mesh of cells can only be applied in 3D because in 2D the extracellular space would become disconnected.

The splitting algorithm requires the solution of a 3D Laplace equation for the entire extracellular domain with Neumann boundary conditions on the cell membranes. With a careful treatment of the boundary conditions, the finite difference discretization produces a symmetric and positive-definite matrix for the linear system associated with the extracellular domain. These properties of the matrix permit the use of a conjugate gradient (CG) iterative solver with an algebraic multigrid (AMG) preconditioner to solve the linear system. Under optimal conditions, the CG+AMG solver should require a constant number of iterations, irrespective of the number of degrees of freedom in the linear system, resulting in a solution time that grows linearly with the size of the linear system [38].

We used the CG and AMG implementations provided by ViennaCL [39] in our simulations, as ViennaCL supports shared-memory parallelization via OpenMP. The termination criterion | A x − b | / | b | < 10 − 5 was used. The coarsening method inside AMG was set to maximum independent set (MIS) [40], and the interpolation method was set to smoothed aggregation. Additionally, the Jacobi smoother weight was set to 0.85 as this value seemed to result in the lowest number of iterations.

On each mesh node along the membrane, the ODE system given in Step 1 of Algorithm 2 is solved with a forward Euler scheme using a substepping time step Δ t * = 1   μ s . We also implemented a first-order generalized Rush–Larsen (GRL1) scheme, and found this theoretically more stable scheme to take about twice the computing time of forward Euler. Since the forward Euler scheme produced stable solutions for all the numerical experiments reported in this paper, the GRL1 scheme was not used. Within each time step, the ODE model can be solved independently on each membrane node, hence there is ample parallelism. We used the Gotran code generator [41] to generate C code for solving a single step for both the forward Euler and GRL1 schemes. The auto-generated code was then augmented with OpenMP directives enabling shared-memory parallel execution.

When initializing the memory for the ODE-model’s state variables and parameters, we perform a so-called “first touch” to ensure that all CPU cores work on the parts of the memory that are located closest to the cores in terms of connectivity. This step is important for achieving proper scaling on modern servers where the system memory is typically split into multiple NUMA (non-uniform memory access) domains.

Two hardware testbeds were used to measure the performance of our parallel solver that implements the splitting algorithm. The specifications are listed in Table 5. The first testbed (System A) is a dual-socket server housing two AMD EPYC 7601 32-core CPUs, each configured with 8-channel memory operating at 2666 MT/s. This server has 2 TiB of memory in total. The second testbed (System B) is a quad-socket server with four Intel Xeon Platinum 8,168 24-core CPUs, each configured with 6-channel memory operating at 2,666 MT/s. The total memory size is 384 GiB. Notably, System A has four NUMA domains per socket, ² whereas System B only has one NUMA domain per socket. Both systems support 2-way simultaneous multi-threading (SMT), presenting two logical processors to the operating system for each physical CPU core. Our hardware test systems are arguably in the high end of the present-day CPU server market and thus allow us to explore the limits of what can be achieved in a shared-memory environment. Simulations with an even larger number of cells, or a finer resolution inside each cell, will eventually require an enhanced parallel simulator capable of utilizing multiple nodes on a distributed-memory system. Speedup due to increasing the number of nodes used can be expected as long as the communication overhead does not exceed the computation time per node.

Our solver is mainly written in the Python programming language, but the most performance-critical parts are executed in C code which we call from Python via the ctypes module. The resulting multi-threaded parallel solver permits rapid prototyping while also being reasonably performant. In the results presented below, we have set environment variables to impose a “thread binding” that fixes the mapping of OpenMP threads to the logical processors such that the same core will work on the same part of the memory for every time step. All logical processors were assigned one thread each, resulting in 128 threads for System A and 192 threads for System B.

The mesh resolution was set to Δ x = Δ y = Δ z = 2   μ m , and we use N it = M it = 2 in the splitting algorithm. N E is the number of unknowns in the linear system for the extracellular domain. I t CG is the mean number of preconditioned CG iterations required for the extracellular system. The time spent solving each of the three domains (Extracellular, Membrane and Intracellular) is listed as well as the total across all domains. On the right side, the times are divided by the number of cells to demonstrate the asymptotic scaling.

Tables 6 and 7 show the time usage per time step needed for 16 cells up to 65,536 and 16,348 cells for Systems A and B, respectively (The memory size on System B is not large enough to simulate 256 × 256 cells.) More specifically, the tables report the average time usage per time step for time steps 6–10 in the simulation. As the number of cells grows, we can see the solution time per cell converging toward a constant. The number of AMG-preconditioned CG iterations for the extracellular system grows by about 71% while the number of unknowns increases by more than a factor of 2,800. System B outperforms System A for all the cases, which can be attributed to it having 50% higher theoretical memory bandwidth, 50% more cores, and higher clock speed than System A. Additionally, architectural differences between the CPUs could also have played a role.

In Table 8, we consider the case of a 4 μm resolution for System A, and we are in this case able to perform simulations of up to 512 × 256 cells. Note that in order to use a 4 μm resolution, the tissue geometry has been slightly altered. The distance from the cells to the extracellular boundary in the z-direction was extended to 8 μm, Ω O was of size 104 μm × 16 μm × 16 μm, and Ω W , Ω E , Ω S , and Ω N had size 8 μm × 8 μm × 8 μm (see Figure 3).

The mesh resolution was set to Δ x = Δ y = Δ z = 4   μ m , and we use N it = M it = 2 in the splitting algorithm. N E is the number of unknowns in the linear system for the extracellular domain. I t CG is the mean number of preconditioned CG iterations required for the extracellular system. The time spent solving each of the three domains (Extracellular, Membrane and Intracellular) is listed as well as the total across all domains. On the right side, the times are divided by the number of cells to demonstrate the asymptotic scaling.