A cut finite element method for spatially resolved energy metabolism models in complex neuro-cell morphologies with minimal remeshing

Abstract

A thorough understanding of brain metabolism is essential to tackle neurodegenerative diseases. Astrocytes are glial cells which play an important metabolic role by supplying neurons with energy. In addition, astrocytes provide scaffolding and homeostatic functions to neighboring neurons and contribute to the blood–brain barrier. Recent investigations indicate that the complex morphology of astrocytes impacts upon their function and in particular the efficiency with which these cells metabolize nutrients and provide neurons with energy, but a systematic understanding is still elusive. Modelling and simulation represent an effective framework to address this challenge and to deepen our understanding of brain energy metabolism. This requires solving a set of metabolic partial differential equations on complex domains and remains a challenge. In this paper, we propose, test and verify a simple numerical method to solve a simplified model of metabolic pathways in astrocytes. The method can deal with arbitrarily complex cell morphologies and enables the rapid and simple modification of the model equations by users also without a deep knowledge in the numerical methods involved. The results obtained with the new method (CutFEM) are as accurate as the finite element method (FEM) whilst CutFEM disentangles the cell morphology from its discretisation, enabling us to deal with arbitrarily complex morphologies in two and three dimensions.

Introduction

We propose to test and verify a simple numerical framework to solve a simplified model of metabolic pathways, representative of cellular metabolism in the brain. Metabolic models can aid the understanding of cell behaviour. Most metabolic models involve the solution of a system of ODEs [1] leaving open the question of how the geometry and spatial organization affect cell behaviour. This work presents a method that is a first step towards extending existing models to answer this question. The method, based on recent developments in unfitted finite element methods is general, and can deal with an arbitrary number of coupled reaction diffusion equations and handle complex cell morphologies. Thanks to the versatility of the automatic code infrastructure provided by FEniCS, the code is well-suited to newcomers to finite element methods interested in modelling biological systems.

To test the new framework, we address complex cell geometries in two and three dimensions and compare the method to the standard finite element method in terms of accuracy. The results obtained with the new method CutFEM are as accurate as the finite element method (FEM) whilst CutFEM disentangles the cell morphology from its discretisation, enabling to deal with arbitrarily complex morphologies, including kinks and triple junctions, in two and three dimensions.

A thorough understanding of brain metabolism is essential to tackle neurodegenerative diseases [2, 3]. Astrocytes are glial cells which play an important metabolic role by supplying neurons with energy [4]. These cells also provide scaffolding and help repair neighboring neurons, where they maintain balanced ionic concentrations (homeostasis) and contribute to the blood–brain barrier by preventing the diffusion of large molecules into the brain.

Recent investigations show that the morphology of astrocytes, which can be complex (see Fig. 1), impacts upon their function, in particular the efficiency with which these cells metabolise nutrients and provide neurons with energy [5]. Modelling and simulation could be effective in furthering our understanding brain metabolism by enabling to test biological hypotheses and investigate the relative importance of model parameters on quantities of interest to biologists.

Fig. 1

The complex morphology of an astrocyte from a human post mortem sample obtained by fluorescent (GFAP) super-resolution light microscopy [63] and reconstructed by the machine learning based tool MicMac [5]

The challenges involved with modelling and simulating metabolic activity in cells include:

• 1. Building a representative model of the metabolic pathways, usually a set of reaction–diffusion equations;

• 2. Identifying the parameters for these partial differential equations (PDEs) as well as the sensitivity of the system to these parameters [6,7,8,9,10,11];

• 3. Discretising the complex and evolving geometries of the cells;

• 4. Solving this set of PDEs on these complex domains;

• 5. Ensuring the accuracy of the solution by measuring discretisation error [12,13,14];

• 6. Ensure the usability of the numerical framework by non-experts.

This last point aims at simplifying multi-disciplinary interactions between computational, data science and domain experts, it is also becoming increasingly important to devise numerical frameworks which can be used and enhanced without being a computational science expert. To do so, open-source frameworks such as FEniCS [15,16,17], getFEM [18], FreeFEM [19], Deal.ii [20] are all possible candidates. These open-source frameworks enable the user to write models in a language which is natural to them and requires minimal interaction with technical details associated with well-established numerical methods.

Concerning metabolic activities, some models focused their attention to the metabolic pathways which they describe as a series of chemical reactions that enable the synthesis and breakdown of molecules as [21,22,23,24].

So far, very little attention has been paid to the role of cell morphology on the modelling [25, 26]. However, cell morphology is known to be important to describe the state of the cell. For example, [27] show the importance of cell morphology on the development of Alzheimer’s disease.

In this paper, we tackle points 1., 3., 4., and 6., above, and leave points 2. and 5. for further communications, as well as the extension to evolving domains.

Model We present a simple model of energy metabolism which takes into account the main pathways of the metabolic process in a single cell. Here, we focus on an astrocyte as a specific cell, but the developments are general. In particular, we include glycolysis, Lactate Dehydrogenase, TCA cycle and basal cellular activity explicitly in the cell model. Each pathway is described by a chemical reaction leading to a coupled reaction diffusion system.

Cell geometry Our ultimate goal is to enable the use of microscopy to produce input geometries for our computational framework and include the evolution of the cellular domain. To enable this, we define cell geometries using signed distance functions, also known as level-set functions [28, 29], which were successfully deployed in modeling other biological phenomena with moving interfaces, such as [30, 31].

Discretisation

The solution of the set of coupled reaction diffusion equations has two characteristics, which create challenges for the standard finite element method:

• Local gradients in metabolite concentrations;

• Discontinuities across cell boundaries.

Standard finite element method would require to create a mesh that fits the boundary of the object. Even though much progress has been made in meshing technology [32], the advantages of methods that separate the geometry from the object are very appealing for our ultimate objective of describing a cell evolving in time. Classic FEM would require to build a mesh conforming the object at each time step leading to very high computational drawbacks.

Enriched finite element methods such as the partition of unity FEM (PUFEM) [33, 34], the generalized finite element method (GFEM) [35] and the extended finite element method [36] are ideal to tackle these two challenges [30, 37].

Indeed, these methods enable the local enrichment of the discrete solution space with known features about the solution, including discontinuities and sharp gradients or singularities. This makes it possible to handle arbitrarily complex geometries quasi-independently of the mesh (see Fig. 2).

Fig. 2

a Standard FEM requires the mesh to be conform the boundary of the domain $\partial \Omega$ b CutFEM embeds the implicitly defined interface $\partial \Omega$ into a background mesh, in yellow the set of ghost penalty facets $\mathcal {F}_G$

Nonetheless, without preconditioning or special treatment [38,39,40,41,42,43,44] enriched finite elements cannot natively deal with arbitrarily complex geometries because of particular geometrical limit cases (interfaces passing close to a degree of freedom) and ill-conditioning stemming from linear dependencies due to complex enrichment functions acting upon large parts of the computational domain.

CutFEM [45,46,47,48,49,50,51] is an extension of XFEM that naturally addresses limit cases associated with complex geometries, and lends itself to image-based simulations. Moreover, we use the libCutFEM library, which is a cut finite element extension of the open-source framework of the FEniCS Project [15,16,17, 47]. FEniCS offers a highly flexible and easy way of transforming models expressed as partial differential equations into numerical methods based on the finite element method through the Unified Form Language (UFL) [52], a domain-specific programming language for writing the weak form of partial differential equations. The UFL specification of the finite element problem is then automatically translated by the other components of the FEniCS Project (DOLFIN [16] and the FEniCS Form Compiler (FFC) [15]) into high-performance specific C++ code with little or no user intervention.

Solution scheme To solve the set of coupled, time-dependent and non-linear PDEs, we first discretise in time using a standard implicit time stepping scheme. The non-linear equation is solved using a Newton–Raphson scheme. The Jacobian, required for the Newton–Raphson scheme, is calculated automatically at the symbolic level by FEniCS. This greatly eases the implementation from the user’s perspective, as deriving and implementing the consistent Jacobian manually can be a tedious and error prone task. The resulting linear systems at each step of the Newton–Raphson algorithm is solved using standard linear solvers from the PETSc library [53].

Usability and extensibility of the framework

Thanks to the flexibility of FEniCS, the framework is generic and usable by biologists without a background or specific training in Computational Sciences or Applied Mathematics. The code used to create the two of the examples is freely online accessible (see “Availability of data and materials”) and can be adapted to other problems straightforwardly.

The paper is organized as follows. In “The problem formulation” section, we present the biological model for energy metabolism followed by the governing equation and its weak form. The FEM and CutFEM discretization can be found in “Discretization” section. In “Implementation” section a detailed explanation of how implement our model in CutFEM is presented. Our numerical results are introduced in “Numerical results” section and the conclusions follow in “Conclusion and discussion” section.

The problem formulation

The aim of this section is to introduce the energy metabolism model, which can be expressed as a set of partial differential equations, in its strong and weak form.

Basic model for energy metabolism

The scope of this model is to isolate conceptually the essential mathematical properties of the metabolic processes of a cell, the mechanism that generate energy for cellular activities from the synthesis and breakdown of nutrients [1].

The simplified model of metabolic pathway in a cell is sketched in Fig. 3 and described by the following reactions (see “Abbreviations” for details)

$$\begin{aligned} \mathrm {HXK}&:= \mathrm {GLC}+ 2\,\mathrm {ATP}\rightarrow 2\,\mathrm {ADP}+ 2 \mathrm {GLY} \end{aligned}$$

(1)

$$\begin{aligned} \mathrm {PYRK}&:= \mathrm {GLY}+2\,\mathrm {ADP}\rightarrow 2\,\mathrm {ATP}+ \mathrm {PYR}\end{aligned}$$

(2)

$$\begin{aligned} \mathrm {LDH}&:= \mathrm {PYR}\rightarrow \mathrm {LAC}\end{aligned}$$

(3)

$$\begin{aligned} \mathrm {Mito}&:= \mathrm {PYR}\rightarrow 28 \,\mathrm {ATP}\end{aligned}$$

(4)

$$\begin{aligned} \mathrm {act}&:= \mathrm {ATP}\rightarrow \mathrm {ADP}. \end{aligned}$$

(5)

Fig. 3

The main pathways we consider in our Simplified Model of Energy Metabolism of a cell. The glucose ($\mathrm {GLC}$) enters inside the cell into the cytosol and takes part into the glycolysis composed by the two reactions $\mathrm {HXK}$ and $\mathrm {PYRK}$. The products of the glycolysis, $\mathrm {PYR}$ and $\mathrm {ATP}$, are then used into $\mathrm {LDH}$ reaction, that generates $\mathrm {LAC}$, into $\mathrm {act}$ reaction, that describes the cellular activity producing $\mathrm {ADP}$, and inside the mitochondria where the reaction denoted $\mathrm {Mito}$ happens generating the major source of $\mathrm {ATP}$ for the cell

The pathway starts when the molecules derived from food (nutrients) enter the cytosol of the cell, and the process called glycolysis starts the breakdown of glucose ($\mathrm {GLC}$) into two molecules of pyruvate ($\mathrm {PYR}$). The glycolysis process can be split, to our modeling purpose, into two main chemical reactions: ATP-consuming (1) and ATP-producing (2). The pyruvate produced by glycolysis can, then, be converted into lactate ($\mathrm {LAC}$) by the enzyme lactate dehydrogenase ($\mathrm {LDH}$) simplified with the chemical reaction (3) or enter mitochondria and used to produce $\mathrm {ATP}$ through the TCA cycle shown in reaction (4).

Last, we take into account the activity of the cell which uses $\mathrm {ATP}$ for its own sustenance represented by the chemical reaction (5).

In order to facilitate our model we consider that the backward reactions are negligible. Moreover, we consider that the enzymes that catalyzed the chemical reactions are located in a specific region of the cellular domain.

Strong formulation of governing equations

Let $\Omega$ be an open and bounded subset of $\mathbb {R}^d$ (d = 2 or 3) with Lipschitz-continous boundary and we denote the concentration of the chemical species using the bracket notation, $[\,\cdot \,]$, as a function

$$\begin{aligned} {[}\, \cdot \,]:\Omega \times [0,T] \rightarrow \mathbb {R}. \end{aligned}$$

Mathematically, we can express the sequence of reactions Eqs. (1)–(5) with a coupled system of semi-linear parabolic reaction diffusion equations [54]. We consider that the species involved are subject to diffusion in the domain and we denote $D_{[\, \cdot \,]} > 0$ the diffusive constants. The reactions obey to the law of mass action [55] where $K_{\mathrm {HXK}},K_{\mathrm {PYRK}},K_{\mathrm {LDH}}, K_{\mathrm {Mito}}$ and $K_{\mathrm {act}}$ are the rate constants and we introduce Gaussian functions, indicated with $\mathcal {G}_{\mathrm {HXK}}(x_0,\sigma )$, $\mathcal {G}_{\mathrm {PYRK}}(x_0,\sigma )$, $\mathcal {G}_{\mathrm {LDH}}(x_0,\sigma )$, $\mathcal {G}_{\mathrm {Mito}}(x_0,\sigma )$ and $\mathcal {G}_{\mathrm {act}}(x_0,\sigma )$, to locate the region where the reactions are happening

$$ \mathcal {G}_\cdot = \mathcal {G}(x_{0,\cdot },\sigma _\cdot ) =\frac{1}{\sqrt{2\pi \sigma _\cdot ^2}}\exp {-\frac{(x-x_{0,\cdot })^2}{2\sigma _\cdot ^2}},$$

where $x_{0,\cdot }\in \Omega _\cdot$ and $\sigma \in \mathbb {R}$. Eventually, to represent the entrance of the glucose inside the cytosol we define a source term function $f: A \times [0,T]\rightarrow \mathbb {R}$, where A is a subset of $\Omega$ as

$$\begin{aligned} f(x,t)={\left\{ \begin{array}{ll} \alpha \in \mathbb {R} \quad {\text {if} \quad (x,t)\in A\times [0,1],}\\ 0 \quad \text {otherwise}. \end{array}\right. } \end{aligned}$$

The strong form of the reaction diffusion system is then, finding the concentrations $[\mathrm {GLC}](x,t)$, $[\mathrm {ADP}](x,t)$, $[\mathrm {ATP}](x,t)$, $[\mathrm {GLY}](x,t)$, $[\mathrm {PYR}](x,t)$ and $[\mathrm {LAC}](x,t)$ for all $x \in \Omega$ and $t\in [0,T]$ such that

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial [\mathrm {GLC}]}{\partial t}=&{}D_{[\mathrm {GLC}]}\nabla ^2[\mathrm {GLC}] - K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}+ f\\ \frac{\partial [\mathrm {ATP}]}{\partial t}=&{}D_{[\mathrm {ATP}]}\nabla ^2[\mathrm {ATP}] -2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}+2K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}} + 28 K_{\mathrm {Mito}}[\mathrm {PYR}]\mathcal {G}_{\mathrm {Mito}}\\ &{}- K_{\mathrm {act}}[\mathrm {ATP}]\mathcal {G}_{\mathrm {act}}\\ \frac{\partial [\mathrm {ADP}]}{\partial t}=&{}D_{[\mathrm {ADP}]}\nabla ^2[\mathrm {ADP}] +2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}-2K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}} +K_{\mathrm {act}}[\mathrm {ATP}]\mathcal {G}_{\mathrm {act}}\\ \frac{\partial [\mathrm {GLY}]}{\partial t}=&{}D_{[\mathrm {GLY}]}\nabla ^2[\mathrm {GLY}] +2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}-K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\\ \frac{\partial [\mathrm {PYR}]}{\partial t}=&{}D_{[\mathrm {PYR}]}\nabla ^2[\mathrm {PYR}] +K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\\ &{}-K_{\mathrm {LDH}} [\mathrm {PYR}] \mathcal {G}_{\mathrm {LDH}} - K_{\mathrm {Mito}}[\mathrm {PYR}]\mathcal {G}_{\mathrm {Mito}}\\ \frac{\partial [\mathrm {LAC}]}{\partial t}=&{}D_{[\mathrm {LAC}]}\nabla ^2[\mathrm {LAC}] + K_{\mathrm {LDH}} [\mathrm {PYR}] \mathcal {G}_{\mathrm {LDH}}.\\ \end{array}\right. } \end{aligned}$$

(6)

The system is completed with homogeneous Neumann boundary conditions on $\partial \Omega$ and initial conditions at time $t=0$

$$\begin{aligned} {\left\{ \begin{array}{ll} {[}\mathrm {GLC}](x,t=0)=a_0{(x)} \quad x\in \Omega \\ {[}\mathrm {ATP}](x,t=0)=b_0{(x)} \quad x\in \Omega \\ {[}\mathrm {ADP}](x,t=0)=c_0{(x)} \quad x\in \Omega \\ {[}\mathrm {GLY}](x,t=0)=d_0{(x)} \quad x\in \Omega \\ {[}\mathrm {PYR}](x,t=0)=e_0{(x)} \quad x\in \Omega \\ {[}\mathrm {LAC}](x,t=0)=f_0{(x)} \quad x\in \Omega . \end{array}\right. } \end{aligned}$$

(7)

where $a_0,b_0,c_0,d_0,e_0,f_0: \Omega \rightarrow \mathbb {R}$.

Precise questions on global existence of solutions of reaction–diffusion systems are still an open problem, we refer the reader to [56]. In this work, we consider all diffusion coefficients $D_{[\cdot ]}$ equals, ensuring that a global solution of the system (6) with initial condition (7) exists.

Weak formulation of governing equations

In this section we convert the strong form of the PDEs in Eq. (6) into a corresponding weak form. This is necessary step in order to discretise both the FEM and CutFEM methods. For further details see e.g. [57].

We define the standard Hilbert space $V={v\in H^1(\Omega )}$ and we denote with $\overline{V}= V \times V \times V \times V \times V \times V$ the product space.

The weak form of system (6) can be found multiplying each equation by a test function $v_1,v_2,v_3,v_4,v_5,v_6\in \overline{V}$ and integrating over the space domain $\Omega$

$$\begin{aligned} {\left\{ \begin{array}{ll} \int _\Omega \frac{\partial [\mathrm {GLC}]}{\partial t}\, v_1 \, \text {d}x =&{} \int _\Omega D_{[\mathrm {GLC}]}\nabla ^2[\mathrm {GLC}]\, v_1 \, \text {d}x - \int _\Omega K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\, v_1 \, \text {d}x\\ &{}+ \int _\Omega f\, v_1\, \text {d}x\\ \int _\Omega \frac{\partial [\mathrm {ATP}]}{\partial t}\, v_2\, \text {d}x=&{} \int _\Omega D_{[\mathrm {ATP}]}\nabla ^2[\mathrm {ATP}]\, v_2\, \text {d}x -\int _\Omega 2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\, v_2 \, \text {d}x\\ &{}+ \int _\Omega 2K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\, v_2 \, \text {d}x\\ &{}+ \int _\Omega 28 K_{\mathrm {Mito}}[\mathrm {PYR}]\mathcal {G}_{\mathrm {Mito}} v_2 \, \text {d}x - \int _\Omega K_{\mathrm {act}}[\mathrm {ATP}]\mathcal {G}_{\mathrm {act}} v_2\, \text {d}x\\ \int _\Omega \frac{\partial [\mathrm {ADP}]}{\partial t}\, v_3\, \text {d}x=&{} \int _\Omega D_{[\mathrm {ADP}]}\nabla ^2[\mathrm {ADP}]\, v_3\, \text {d}x +\int _\Omega 2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\, v_3 \, \text {d}x\\ &{}- \int _\Omega 2K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\, v_3 \, \text {d}x+ \int _\Omega K_{\mathrm {act}}[\mathrm {ATP}]\mathcal {G}_{\mathrm {act}}\, v_3 \, \text {d}x\\ \int _\Omega \frac{\partial [\mathrm {GLY}]}{\partial t}\, v_4\, \text {d}x =&{}\int _\Omega D_{[\mathrm {GLY}]}\nabla ^2[\mathrm {GLY}]\, v_4\, \text {d}x +\int _\Omega 2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\, v_4 \, \text {d}x\\ &{}- \int _\Omega K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\, v_4\, \text {d}x\\ \int _\Omega \frac{\partial [\mathrm {PYR}]}{\partial t}\, v_5\, \text {d}x=&{}\int _\Omega D_{[\mathrm {PYR}]}\nabla ^2[\mathrm {PYR}]\, v_5\, \text {d}x +\int _\Omega K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\, v_5 \, \text {d}x\\ &{}-\int _\Omega K_{\mathrm {LDH}} [\mathrm {PYR}] \mathcal {G}_{\mathrm {LDH}}\, v_5 \, \text {d}x - \int _\Omega K_{\mathrm {Mito}}[\mathrm {PYR}]\mathcal {G}_{\mathrm {Mito}}\, v_5\, \text {d}x\\ \int _\Omega \frac{\partial [\mathrm {LAC}]}{\partial t}\, v_6\, \text {d}x=&{}\int _\Omega D_{[\mathrm {LAC}]}\nabla ^2[\mathrm {LAC}]\, v_6\, \text {d}x+\int _\Omega K_{\mathrm {LDH}} [\mathrm {PYR}] \mathcal {G}_{\mathrm {LDH}}\, v_6\, \text {d}x.\\ \end{array}\right. } \end{aligned}$$

(8)

Since $D_{[\cdot ]}$ is constant, we apply integration by parts to the second-order spatial derivative

$$\begin{aligned} -\int _\Omega D_{[\, \cdot \,]} \nabla ^2 [\, \cdot \,]\, v_i\, \text {d}x = \int _\Omega D_{[\, \cdot \,]} \nabla [\, \cdot \,]\, \nabla v_i\, \text {d}x -\int _{\partial \Omega } D_{[\, \cdot \,]} \frac{\partial [\, \cdot \,]}{\partial n} v_i \text {d}s. \end{aligned}$$

(9)

Since we have specified pure Neumann boundary conditions for each concentration species, the boundary terms vanish. Leading to

$$\begin{aligned} -\int _\Omega D_{[\, \cdot \,]} \nabla ^2 [\, \cdot \,]\, v_i\, \text {d}x = \int _\Omega D_{[\, \cdot \,]} \nabla [\, \cdot \,]\, \nabla v_i\, \text {d}x. \end{aligned}$$

(10)

We can, then, substitute Eq. (10) into Eq. (8) and use the following compact notation

$$\begin{aligned} (u,v) = \int _\Omega uv \, \text {d}x. \end{aligned}$$

(11)

allowing us to state the weak form of the problem problem as: we seek the solutions $[\mathrm {GLC}]$,$[\mathrm {ATP}]$, $[\mathrm {ADP}]$,$[\mathrm {GLY}]$, $[\mathrm {PYR}]$, $[\mathrm {LAC}]$ in the space $\overline{V}$ such that for all $t\in [0,T]$ and for all test functions $v_1,v_2,v_3,v_4,v_5,v_6\in \overline{V}$

$$\begin{aligned} {\left\{ \begin{array}{ll} \Big (\frac{\partial [\mathrm {GLC}]}{\partial t},v_1\Big )=&{}\Big (-D_{[\mathrm {GLC}]}\nabla [\mathrm {GLC}],\nabla v_1\Big )+\Big ( - K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}+ f,v_1\Big )\\ \Big (\frac{\partial [\mathrm {ATP}]}{\partial t},v_2\Big )=&{}\Big (-D_{[\mathrm {ATP}]}\nabla [\mathrm {ATP}],\nabla v_2\Big ) +\Big (-2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}+2K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}} + 28 K_{\mathrm {Mito}}[\mathrm {PYR}]\mathcal {G}_{\mathrm {Mito}}\\ &{}- K_{\mathrm {act}}[\mathrm {ATP}]\mathcal {G}_{\mathrm {act}},v_2\Big )\\ \Big (\frac{\partial [\mathrm {ADP}]}{\partial t},v_3\Big )=&{}\Big (-D_{[\mathrm {ADP}]}\nabla [\mathrm {ADP}],\nabla v_3\Big ) +\Big (+2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}-2K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}+K_{\mathrm {act}}[\mathrm {ATP}]\mathcal {G}_{\mathrm {act}},v_3\Big )\\ \Big (\frac{\partial [\mathrm {GLY}]}{\partial t},v_4\Big )=&{}\Big (-D_{[\mathrm {GLY}]}\nabla [\mathrm {GLY}],\nabla v_4\Big ) +\Big ( +2K_{\mathrm {HXK}}[\mathrm {GLC}][\mathrm {ATP}]^2\mathcal {G}_{\mathrm {HXK}}\\ &{}-K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}},v_4\Big )\\ \Big (\frac{\partial [\mathrm {PYR}]}{\partial t},v_5\Big )=&{}\Big (-D_{[\mathrm {PYR}]}\nabla [\mathrm {PYR}],\nabla v_5\Big ) +\Big (+K_{\mathrm {PYRK}}[\mathrm {ADP}]^2[\mathrm {GLY}]\mathcal {G}_{\mathrm {PYRK}}\\ &{}-K_{\mathrm {LDH}} [\mathrm {PYR}] \mathcal {G}_{\mathrm {LDH}} - K_{\mathrm {Mito}}[\mathrm {PYR}]\mathcal {G}_{\mathrm {Mito}},v_5\Big )\\ \Big (\frac{\partial [\mathrm {LAC}]}{\partial t},v_6\Big )=&{}\Big (-D_{[\mathrm {LAC}]}\nabla [\mathrm {LAC}],\nabla v_6\Big )+\Big ( + K_{\mathrm {LDH}} [\mathrm {PYR}] \mathcal {G}_{\mathrm {LDH}},v_6\Big ).\\ \end{array}\right. } \end{aligned}$$

(12)

In the subsequent section we discuss the discretization of Eq. (12) using the standard finite element method and the CutFEM.

Discretization

In order to solve Eq. (12), we must discretise it in both space and time. Discretisation is a process by which continuous mathematical objects (e.g. $v_1$) are transformed into a discrete counterpart that can be manipulated on a computer. We choose to discretise in space using the classical Finite Element Method [58] (FEM) and then using the Cut Finite Element Method (CutFEM) [47]. The FEM results serve as a baseline for comparison of the CutFEM method. For both FEM and CutFEM we discretize in time using a standard finite difference method.

The important distinction between FEM and CutFEM from the point of view of the user is that the FEM requires an extra mesh generation step; a mesh must be generated that conforms to the boundary of the domain, before the simulation can take place. This can be a difficult task, as the mesh must be of sufficiently good quality to ensure an accurate solution, while still conforming to the boundary. In contrast, CutFEM removes the need for a conforming mesh generation step. The boundary of the domain is described implicitly as a level set function that can be extracted directly from e.g. processed image data, or using constructive solid geometry (CSG) [59]. The promise of CutFEM is that discretization of geometry can be performed automatically without a mesh generation step that often requires lengthy manual intervention.

We recognise that both the Finite Element Method, and even more so, the CutFEM are highly technical and take some time to understand. The point of this section then is not to give a full and detailed exposition of both of these methods. Instead we aim show a precise derivation of the discrete weak forms and then in the subsequent section we show their translation into the FEniCS Project domain specific language, called the Unified Form Language (UFL) [52]. In practice, if the user can convert their problem into a discrete weak formulation then the the FEniCS Project and the libCutFEM library can automatically perform all of the subsequent discretization steps. For full details of how this process takes places the reader is referred to [16, 47].

FEM

To discretise our problem in time we choose a finite difference approximation, specifically, the backward Euler method [60]. First, we discretize the time interval [0, T], denoting $\Delta t$ the size of the timestep. We denote with subscript n a concentration of (possibly multiple) chemical species at time $t_n$, where $0 \le n < N$ with $N = T/\Delta t$ is an integer that counts the time steps. The backward Euler method then approximates the continuous time derivative as

$$\begin{aligned} \frac{\partial }{\partial t} [\, \cdot \,] \approx \frac{[\,\cdot \,]_{n+1} - [\,\cdot \,]_n}{\Delta t}. \end{aligned}$$

(13)

We discretize the spatial domain $\Omega$ using a triangulation $\mathcal {T}_h$ that conforms (matches) the exact boundary $\partial \Omega$. On this mesh we define the space of piece wise Lagrange finite elements of degree one as V and with $\overline{V}$ the product $\overline{V} = V \times V \times V \times V \times V \times V$. The space $\overline{V}$ can then be used to discretize the weak form equations (12).

We denote $U=([\mathrm {GLC}],[\mathrm {ATP}],[\mathrm {ADP}],[\mathrm {GLY}],[\mathrm {PYR}],[\mathrm {LAC}])$ the vector of all the concentration solutions and $v=(v_1,v_2,v_3,v_4,v_5,v_6) \in \overline{V}$ the corresponding vector of test functions. The initial conditions are denoted $[U]_0$. We then solve a sequence of problems: find $[U]_{n + 1} \in \overline{V}$ for $n = 0, \ldots , N - 1$ such that

$$\begin{aligned} f_h([U]_{n + 1}, [U]_{n} ;v) = 0 \quad \forall v\in \overline{V}, \end{aligned}$$

(14)

where

$$\begin{aligned} f_h([U]_{n+1}, [U]_{n} ;v)&= \int _\Omega \Big (\Delta t^{-1} ([\mathrm {GLC}]_{n+1}-[\mathrm {GLC}]_n)v_1 + D_{[\mathrm {GLC}]}\nabla [\mathrm {GLC}]_{n+1}\nabla v_1 +\nonumber \\&-fv_1 +K_{\mathrm {HXK}}[\mathrm {GLC}]_{n+1}[\mathrm {ATP}]_{n+1}^2 v_1\mathcal {G}_{\mathrm {HXK}}+\nonumber \\&+\Delta t^{-1} ([\mathrm {ATP}]_{n+1}-[\mathrm {ATP}]_n)v_2 + D_{[\mathrm {ATP}]}\nabla [\mathrm {ATP}]_{n+1}\nabla v_2 +\nonumber \\&+2K_{\mathrm {HXK}}[\mathrm {GLC}]_{n+1}[\mathrm {ATP}]_{n+1}^2 v_2\mathcal {G}_{\mathrm {HXK}}+\nonumber \\&-2K_{\mathrm {PYRK}}[\mathrm {GLY}]_{n+1}[\mathrm {ADP}]_{n+1}^2 v_2\mathcal {G}_{\mathrm {PYRK}}+\nonumber \\&-28 K_{\mathrm {Mito}}[\mathrm {PYR}]_{n+1} v_2\mathcal {G}_{\mathrm {Mito}}+ K_{\mathrm {act}}[\mathrm {ATP}]_{n+1} v_2\mathcal {G}_{\mathrm {act}}+\nonumber \\&+\Delta t^{-1} ([\mathrm {ADP}]_{n+1}-[\mathrm {ADP}]_n)v_3 + D_{[\mathrm {ADP}]}\nabla [\mathrm {ADP}]_{n+1}\nabla v_3\nonumber \\&-2K_{\mathrm {HXK}}[\mathrm {GLC}]_{n+1}[\mathrm {ATP}]_{n+1}^2 v_3\mathcal {G}_{\mathrm {HXK}} +\nonumber \\&+2K_{\mathrm {PYRK}}[\mathrm {GLY}]_{n+1}[\mathrm {ADP}]_{n+1}^2 v_3\mathcal {G}_{\mathrm {PYRK}}- K_{\mathrm {act}}[\mathrm {ATP}]_{n+1} v_3\mathcal {G}_{\mathrm {act}}\nonumber \\&+\Delta t^{-1} ([\mathrm {GLY}]_{n+1}-[\mathrm {GLY}]_n)v_4 + D_{[\mathrm {GLY}]}\nabla [\mathrm {GLY}]_{n+1}\nabla v_4 +\nonumber \\&-2K_{\mathrm {HXK}}[\mathrm {GLC}]_{n+1}[\mathrm {ATP}]_{n+1}^2 v_4\mathcal {G}_{\mathrm {HXK}} +\nonumber \\&+K_{\mathrm {PYRK}}[\mathrm {GLY}]_{n+1}[\mathrm {ADP}]_{n+1}^2 v_4\mathcal {G}_{\mathrm {PYRK}}\nonumber \\&+\Delta t^{-1} ([\mathrm {PYR}]_{n+1}-[\mathrm {PYR}]_n)v_5 + D_{[\mathrm {PYR}]}\nabla [\mathrm {PYR}]_{n+1}\nabla v_5 +\nonumber \\&- K_{\mathrm {PYRK}}[\mathrm {GLY}]_{n+1}[\mathrm {ADP}]_{n+1}^2 v_5\mathcal {G}_{\mathrm {PYRK}} + \nonumber \\&+K_{\mathrm {LDH}}[\mathrm {PYR}]_{n+1} v_5\mathcal {G}_{\mathrm {LDH}} + K_{\mathrm {Mito}}[\mathrm {PYR}]_{n+1} v_5\mathcal {G}_{\mathrm {Mito}}\nonumber \\&+\Delta t^{-1} ([\mathrm {LAC}]_{n+1}-[\mathrm {LAC}]_n)v_6 + D_{[\mathrm {LAC}]}\nabla [\mathrm {LAC}]_{n+1}\nabla v_6 +\nonumber \\&-K_{\mathrm {LDH}}[\mathrm {PYR}]_{n+1} v_6\mathcal {G}_{\mathrm {LDH}} \Big )\, \text {d} x. \end{aligned}$$

(15)

Since Eq. (14) is a non-linear function of the unknown solution $[U]_{n + 1}$ we choose to use a Newton–Raphson type algorithm to solve it. This Newton–Raphson algorithm requires the computation of the derivative of $f_h$ with respect to U (commonly called the Jacobian of $f_h$). We do not perform this step manually, but instead use the automatic differentiation capabilities of UFL [52], as shown in “Implementation” section.

CutFEM

Instead of discretizing the spatial domain $\Omega$ using a mesh that conforms to the boundary, in CutFEM, the problem domain $\Omega$ is described by a level-set function. The level set function is a scalar function on $\mathbb {R}^d$, such that $\phi (x)<0$ for $x \in \Omega$, $\phi (x)>0$ for $x \notin \Omega$ and $\phi (x)=0$ for $x \in \partial \Omega$. We then cover the domain $\Omega$ by a regular background mesh $\Lambda$ ($\Omega \subseteq \Lambda$) of simple shape, e.g. a box containing $\Omega$ meshed with a uniform triangulation. Let K denote a triangle/tetrahedron in this triangulation. Now, let $\tilde{\mathcal {T}}_h$ be the fictitious domain mesh composed by all elements $K \in \Lambda$ such that $K \cap \Omega \ne 0$ ($\Omega \subseteq \tilde{\mathcal {T}}_h \subseteq \Lambda$).

Furthermore, the union of all elements in $\tilde{\mathcal {T}}_h$ is called the fictitious domain $\tilde{\Omega }$. We denote with

$$\begin{aligned} \mathcal {G}^*_h:=\{K\in \tilde{\mathcal {T}}_h: K\cap \partial \Omega \ne \emptyset \} \end{aligned}$$

the set of elements intersected by the interface, and we define the set of so-called ghost penalty facets [46] (see Fig. 2)

$$\begin{aligned} \mathcal {F}_G=\{F \, \text {facet in } \tilde{\mathcal {T}}_h: F=K\cap K' \quad \text {where} \quad K\in \mathcal {G}^*_h \quad \text {or}\quad K'\in \mathcal {G}^*_h \}. \end{aligned}$$

The stabilisation term introduced in the next section will be applied to this subset of facets. We consider the space

$$\begin{aligned} W=\{ v\in C^0(\tilde{\Omega }): v|_K\in P^1(K), \forall K\in \tilde{\mathcal {T}}_h\} \end{aligned}$$

and the jump gradient is defined for all facet F and $v \in W$ by $[\![ \partial _{n_F}v]\!]=n_F\cdot \nabla v|_K-n_F\cdot \nabla |_{K'}$ where $n_F$ denotes the unit normal to F in fixed but arbitrary direction. We use the same notation as the previous section, denoting the solution $U=([\mathrm {GLC}],[\mathrm {ATP}],[\mathrm {ADP}],[\mathrm {GLY}],[\mathrm {PYR}],[\mathrm {LAC}])$, $v=(v_1,v_2,v_3,v_4,v_5,v_6)$ and the product space $\overline{W}= W \times W \times W \times W \times W \times W$, We then solve a sequence of problems: find $[U]_{n + 1} \in \overline{V}$ for $n = 0, \ldots , n - 1$ such that

$$\begin{aligned} f_h([U]_{n + 1}, [U]_{n};v) + j([U]_{n + 1}, v) = 0 \quad \forall v\in \overline{W} \end{aligned}$$

(16)

where $f_h([U]_{n + 1}, [U]_{n};v)$ is identical to the standard FEM Equation (15). Here, $j([U]_{n + 1}, v)$ denotes the following stabilization terms

$$\begin{aligned} j([U]_{n+1},v)= \sum _{F\in \mathcal {F}_G}&\big ((\gamma h D_{[\mathrm {GLC}]} [\![ \partial _{n_F}[\mathrm {GLC}]_{n+1}]\!], [\![ \partial _{n_F}v_1]\!])_F \end{aligned}$$

(17)

$$\begin{aligned} +&(\gamma h D_{[\mathrm {ATP}]}[\![ \partial _{n_F}[\mathrm {ATP}]_{n+1}]\!], [\![ \partial _{n_F}v_2]\!])_F\end{aligned}$$

(18)

$$\begin{aligned} +&(\gamma h D_{[\mathrm {ADP}]}[\![ \partial _{n_F}[\mathrm {ADP}]_{n+1}]\!], [\![ \partial _{n_F}v_3]\!])_F\end{aligned}$$

(19)

$$\begin{aligned} +&(\gamma h D_{[\mathrm {GLY}]}[\![ \partial _{n_F}[\mathrm {GLY}]_{n+1}]\!], [\![ \partial _{n_F}v_4]\!])_F\end{aligned}$$

(20)

$$\begin{aligned} +&(\gamma h D_{[\mathrm {PYR}]}[\![ \partial _{n_F}[\mathrm {PYR}]_{n+1}]\!], [\![ \partial _{n_F}v_5]\!])_F\end{aligned}$$

(21)

$$\begin{aligned} +&(\gamma h D_{[\mathrm {LAC}]}[\![ \partial _{n_F}[\mathrm {LAC}]_{n+1}]\!], [\![ \partial _{n_F}v_6]\!])_F\big ) \end{aligned}$$

(22)

inspired by [47]. Here, $\gamma$ is a positive penalty parameter. These stabilization terms extend the solution from the physical domain $\Omega$ onto the fictitious domain $\tilde{\Omega }$, it is consistent with the continuous system. They prevent ill-conditioning of the system matrices in case only small parts of $\Omega$ are contained in an element near the boundary $\partial \Omega$. This stabilisation is critical for the robustness and reliability of CutFEM.

Implementation

The finite element method discretization has been implemented using Python with the open source finite element solver DOLFIN from the FEniCS Project, see [15, 16]. The CutFEM discretization has been implemented using Python using the libCutFEM library [47] which builds on top of DOLFIN and the rest of the FEniCS Project. In this section we show parts of the code for the libCutFEM example that highlight the close link between the mathematics and the concrete computer implementation. The standard FEniCS Project implementation is similar so we have chosen to show only the libCutFEM implementation for reasons of brevity. The reader should refer to the free online repository for two full working examples that are around 250 lines of code each. The precise problem setup and results from this example are shown in “Two-dimensional non-Lipschitz domain” section.

We import the dolfin and cutfem Python modules. These two modules contain all of the functionality we need to solve the problem using the CutFEM approach.

We create the background mesh $\Lambda$ and then define a level set function describing the heart-shaped domain.

In the next part of the code, we use two special libCutFEM methods to 1. create a special cut mesh, and 2. a mesh_cutter that intersects the mesh with the level set and computes the distinct sets of cells (i.e. those that are on the inside of the level set and those that are outside).

With these special objects in hand we can define special CutFEM-specific UFL Measure objects that will subsequently allow us to write the residual and the stabilisation weak forms. Simply put, a measure defines regions of the problem mesh (cells, edges, parts inside the level set, and parts outside etc.) on which FEniCS will integrate different parts of a weak form.

We create a finite element function space $\bar{V}$ that we can use to further define the UFL algebraic objects that we need.

Now we have everything that we need to write the weak residual $f_h([U]_{n+1},[U]_n; v)$ in Eq. (14) using the UFL. We use the integration measure dxc which indicates the integration on cells inside the domain $\Omega$ and in parts of the cut cells inside $\Omega$. This form would look identical in the standard FEniCS Project code except that we would use the dx measure that denotes integration over all cells of the mesh.

The stabilization term $j([U]_{n+1},v),$ Eq. (17) is written in UFL using the integration measure dS(1) which integrates on the ghost penalty facets $\mathcal {F}_G$ (Fig. 4).

Fig. 4

We can see the fictitious domain $\tilde{\Omega }$ for the heart shaped function, and the facets marked. In red the set of facets $\mathcal {F}_G$ where we apply the stabilization term $j([U]_{n+1},v)$

For both forms F and j we remark how similar the UFL notation is to the mathematical notation in Eqs. (14) and (17). Calculating the UFL expression for the Jacobian $\mathbb {J}$ can be performed automatically using the derivative function.

The last step before the solver is to use the Composite framework of CutFEM to define the problem on different parts of the mesh domain.

Eventually, we can solve the non-linear and time-dependent problem with the following two nested loops, the outer one for the timestepping, and the inner one for the Newton–Raphson algorithm. The FEniCS Project infrastructure (UFL + FFC + DOLFIN) automatically generates and executes low-level code to assemble the sparse matrices A and b. This linear system is then solved using PETSc and MUMPS.

In the previous piece of code, we used the increment as stopping criteria for the Newton–Raphson algorithm. We refer the reader to the free online repository with an example where the stopping criteria is the tolerance to the residual. In the full example in the free online repository the solution at each timestep is outputted to a VTK file that can be opened with Paraview (https://paraview.org) for visualisation.

Numerical results

In this section we evaluate the accuracy of the CutFEM discretisation scheme for different geometries by comparing the CutFEM solution to the standard FEM solution. In addition, we want to confirm that at the steady state solutions (for large values of time t) predicted by CutFEM tend towards the asymptotic solutions of the associated ordinary differential equation (ODE) system. Then, we investigate the accuracy of CutFEM in comparison with the standard FEM for a simple circular geometry. We increase the complexity of the geometry and we show the ability of CutFEM to straightforwardly solve a test case within a non-Lipschitz domain. Lastly, we consider a three dimensional domain.

The numerical scheme was implemented using the CutFEM library [47] based on FEniCS Project [15,16,17]. The linear systems arising in the numerical experiments are solved using a direct (MUMPS) solver for the two-dimensional examples and a iterative (CG) solver with algebraic multigrid preconditioning (hypre) for the three-dimensional example.

Asymptotic solution ODEs

The aim of this section is to validate our CutFEM implementation highlight that the solution of the reaction–diffusion system tends to the solution of the ODE system associated to the chemical reactions Eqs. (1)–(5) for time going to infinity. The solution of the ODE system is computed in two ways 1. we use the solve $\_$ ivp of the package scipy in python 2. we manually compute the asymptotic solutions, which can be found in Appendix A together with the ODE system in Eq. (24).

For the PDEs, we solve the reaction diffusion system in a circular domain defined in CutFEM using the level set function $\phi (x,y) = (x-4)^2+(y-3)^2 - 25$. The Gaussian parameters locating the chemical reactions, that are shown in Fig. 5a, have been set as following: $\mathcal {G}_{\mathrm {HXK}}(x_0=0.5, y_0=2.0, \sigma =0.1)$, $\mathcal {G}_{\mathrm {PYRK}}(x_0=1.1, y_0=1.2, \sigma =0.1)$, $\mathcal {G}_{\mathrm {LDH}}(x_0=4.0, y_0=5.0,\sigma =0.1)$, $\mathcal {G}_{\mathrm {Mito}}(x_0=4.0,y_0=5.0, \sigma =0.1)$ and $\mathcal {G}_{\mathrm {act}}(x_0=6.0,y_0=6.5, \sigma =0.1)$. Note, we have co-localized the reactions (3) and (4), in order to obtain that $\mathrm {PYR}$ is used equally in the two reactions, such that we can set the parameter $\alpha$ of the ODE asymptotic solutions of Eq. (26) equal to 0.5.

Fig. 5

a The locations of the chemical reactions Eqs. (1)–(5) and the influx of $\mathrm {GLC}$ are shown in the circular domain in order to prove the convergence of the PDEs solutions to the ODEs. Note, the chemical reactions $\mathrm {LDH}$ Eq. (3) and $\mathrm {Mito}$ Eq. (4) have been placed in the same location. b The solutions of the PDEs obtained using CutFEM and integrated over the domain (red line) are plot in comparison with the asymptotic solutions of the ODEs for each concentration (blue line) and with the ODEs solved using scipy. The convergence of the PDEs to the ODEs solutions is reached for $t>200$

The influx of $\mathrm {GLC}$ entering the cell domain has been set equal to 100 until time $t=1$ and is entering the cell into a circular subdomain with radius 0.3 and center the origin.

The rate constants $K_{\cdot }$ of the chemical reactions are set equal to 10.0, in order to avoid that one reaction dominates over, and the diffusive parameters $D_{[\,\cdot \,]}$ of each chemical species is 100.0 to accelerate the convergence to the ODEs solutions. The initial amount of concentrations inside the cell are set to zero except for $[\mathrm {ATP}]$ and $[\mathrm {ADP}]$ that are equal to 1. Final time is 1000, as well as the number of time step. The CutFEM penalty parameter $\gamma$ has been set equal to 0.1. The mesh size is set to a maximum diameter of 0.1744. In this experiment, a finer mesh is not required for proving the convergence to the ODEs solutions.

As initial condition, the asymptotic ODEs use the solution of PDEs at time equal to 1 when the source term of glucose stops, which keeps the total amounts of concentration unchanged, as shown in Appendix A.

To use scipy to solve the ODEs we set all the parameters as the PDEs, the influx of $\mathrm {GLC}$ and the initial conditions are obtained from the PDE ones integrating over the domain.

In order to compare the numerical solution of the system (6) solved using CutFEM, we integrate the solutions of the concentrations over the domain $\Omega$ to obtain the average chemical concentration of each species at each time step. We use the following notation

$$\begin{aligned} \tilde{[ \cdot ]} (t) := \int _\Omega [ \cdot ](t)\, \text {d}x \quad \forall t\in [0,T]. \end{aligned}$$

(23)

Fig. 6

a The chosen location of the chemical reactions Eqs. (1)–(5) and of the influx of $\mathrm {GLC}$ are shown to compare FEM and CutFEMresults in a circular domain. b In the y-axis is plot the integral over the domain at each time step using Eq. (23), giving us the average concentration, of CutFEM and FEM in a circular domain and CutFEM solved in a circular domain with perturbed boundary. The average concentration of each species of FEM (blue line) and CutFEM (red line) are equivalent while the CutFEM solutions of the perturbed domain (black line) share the same results for concentrations $\mathrm {ATP}$, $\mathrm {GLY}$, $\mathrm {PYR}$, $\mathrm {LAC}$ and differs but show the same trend for $\mathrm {GLC}$ and $\mathrm {ADP}$

In Fig. 5b, we plot the solutions obtained using the Formula (23) for the PDEs, the asymptotic ODE solutions obtained using Eq. (26) and the ODE solutions with scipy. As expected from the asymptotic solutions computed in Appendix A, the average concentrations of $\mathrm {GLC}$, $\mathrm {ATP}$, $\mathrm {GLY}$ and $\mathrm {PYR}$ go to zero whilst the products of the system are $\mathrm {LAC}$ and $\mathrm {ADP}$. We can see that the ODE solutions with scipy converge to the asymptotic solution very quickly, and after time $t=200$ also the PDE solutions tend to the same values.

Example one: two-dimensional circular domain

In this section, we assess the accuracy of the CutFEM solution. First, we consider the cell, $\Omega$, as a circle defined by the level set function $\phi (x,y)= (x-4)^2+(y-3)^2 - 25$ and we show that the results obtained with CutFEM are of the same order as those obtained with FEM, where the circular domain is explicitly meshed using package mshr.

To investigate the behaviour of CutFEM for more complex boundaries we consider an irregular cell geometry defined by a perturbed circle represented by the level set function $\phi (x,y) = (x-4)^2+(y-3)^2 - \cos {(4x)}\cos {(5x)} - \sin {(4y)}\cos {(5y)} - 25.$ We highlight how it is easy for CutFEM to work with a more complex shape just by changing the level set function.

In this experiment the glucose source term and the reaction locations and parameters are set as in section , with the exception of the $\mathrm {Mito}$ reaction that is located at $\mathcal {G}_{\mathrm {Mito}}(x_0=4.0, y_0=7.5, \sigma =0.1)$. The test case is represented schematically in Fig. 6a.

The chemical species are allowed to diffuse inside the cell with diffusive constants $D_{[\,\cdot \,]}=1.0$ for each species. The rate constants of the chemical reaction and the initial chemical concentrations are set as in the previous test case, see Section . Final time is 10 and we set the number of time steps as 300. As before, the penalty parameter of CutFEM is set to 0.1. The choice of the mesh size is shown in Table 1 in accordance with the size of the Gaussian parameters $\sigma$.

Table 1 In the table shows the maximum cell diameter, the number of DOFs, and the computational time for FEM and CutFEM in a circular domain and for CutFEM for the perturbed boundary domain

In Fig. 6b, we investigate the average concentrations inside the three domains at each time step using Eq. (23). The average concentrations for each species computed in the circular domain with FEM and CutFEM show equivalent results. The solutions obtained using the perturbed domain show higher average concentration of $\mathrm {GLC}$ and $\mathrm {ADP}$ than the results in the circular domain. We remind the reader that the system solved inside the perturbed domain is not supposed to produce same results as the circular one, however the dynamics show the same trend as shown in Fig. 6b.

For time equal to 5, we plot the solutions of the concentrations for each chemical species in Fig. 7. We notice how the Gaussian functions drive the chemical reactions in the regions defined in Fig. 6a. For example, we can notice that $\mathrm {ATP}$ is consumed in two regions that locate respectively $\mathrm {HXK}$ and $\mathrm {act}$, while $\mathrm {ADP}$ is produced. The $\mathrm {GLC}$ entered the cell is diffusing in the domain and taking part to the $\mathrm {HXK}$ reaction. In comparison, the results of FEM and CutFEM are almost identical.

Fig. 7

Solutions of the concentrations, from top to bottom: $\mathrm {GLC}$, $\mathrm {ATP}$, $\mathrm {ADP}$, $\mathrm {GLY}$, $\mathrm {PYR}$ and $\mathrm {LAC}$, at time $t=5$ in a circular domain solved with FEM (Left), CutFEM (Center) and CutFEM in a perturbed boundary (Right). The source term of $\mathrm {GLC}$ is no longer active, and the quantity of $\mathrm {GLC}$ is spreading through the domain and participating into the reaction $\mathrm {HXK}$. $\mathrm {ATP}$ and $\mathrm {ADP}$ show that in the region where $\mathrm {HXK}$, $\mathrm {PYRK}$ and $\mathrm {act}$ are located, when one is consumed the other is produced. We can notice production of $\mathrm {GLY}$, $\mathrm {PYR}$ and $\mathrm {LAC}$. Comparing the results in the three domain, they are extremely close to each other and visually not distinguishable

To further analyze the plot shown in Fig. 7, we compare in Fig. 8 FEM and CutFEM solutions line $x=y$ at time 5 (Fig. 8a) and final time 10 (Fig. 8b). We can recognize where the reactions are from the bump in the curve. The plot indicates that the results are equivalent.

Fig. 8

We plot the solutions of the concentrations in a circular domain obtained with FEM and CutFEM along the line $x=y$ at time $t=5$ (a) and final time $T=10$ (b). We can notice in the plot the location where the reactions are happening thanks to the bump in the plot. The only small differences detectable are in $\mathrm {LAC}$ at $t=5$ and in $\mathrm {GLC}$ at $t=10$ where CutFEM is slightly below the FEM curve

Our experiments suggest that CutFEM produces similar solutions to the FEM but the computational time needed to solve CutFEM is five times as high as standard FEM. The number of DOFs is smaller for CutFEM compared to FEM, as shown in Table 1.

Concerning the larger computational time requires for CutFEM, the CutFEM implementation is based on FEniCS 1.5.0 released in 2015, while the FEM ones on FEniCS 2019.1.0 released in 2019 with significant improvements and computational optimizations. Since CutFEM relies in large parts on the FEniCS functions, its computational time would be highly reduced when using a newer version of FEniCS. Breaking down the total computational time, we identified that the CutFEM implementation of assembly is significantly slower than the FEniCS implementation of assembly. However, critically, we have verified that both implementations scale at the expected optimal rate O(N) in the number of cells N (results not shown). The larger constant in the O(N) scaling in the CutFEM implementation points to certain computational kernels not optimised for efficiency.

In conclusion, in this section we showed how CutFEM and FEM gives equivalent results for a circular domain. Moreover, we investigated the behaviour of CutFEM with a more complex domain, the perturbed circle.

Two-dimensional non-Lipschitz domain

In this section, we consider an irregular boundary, i.e. where the normal field along the surface is non-smooth. We choose a heart-shaped domain, which has two singularities (this is known as a non-Lipschitz domain).

The heart-shaped is described by the level set function $\phi (x,y) = (y-\sqrt{|x|})^2 - 1 + x^2$. The description of the problem setting is shown in Fig. 9. The gaussian functions are defined as $\mathcal {G}_{\mathrm {HXK}}(x_0=0.1, y_0=-0.5, \sigma =0.1)$, $\mathcal {G}_{\mathrm {PYRK}}(x_0=0.3, y_0=10.0, \sigma =0.1)$, $\mathcal {G}_{\mathrm {LDH}}(x_0=-0.5, y_0=0.5,\sigma =0.1)$, $\mathcal {G}_{\mathrm {Mito}}(x_0=0.5,y_0=0.7,\sigma =0.1)$ and $\mathcal {G}_{\mathrm {act}}(x_0=0.0,y_0=0.9,\sigma =0.1)$. We consider a source term for $\mathrm {GLC}$ active within a disk of radius 0.3, centered at the cusp and take the $\mathrm {GLC}$ influx to be 100.

Fig. 9

Location of the chemical reactions Eqs. (1)–(5) and the source term of $\mathrm {GLC}$ in the heart shaped domain

The diffusive constants, the reaction rate and the initial concentrations are as in “Example one: two-dimensional circular domain” Section, as well as the penalty parameter for the stabilization term.

We solve the reaction diffusion system (6) using CutFEM, where we have set a background mesh with a maximum cell diameter equal to 0.01803, the number of DOFs is 129,606 and the total time to run the simulation was approximately 4 h.

Figures 10 and 11 show the solutions obtained in the heart-shaped domain for each chemical species at the initial and final time. At the initial time, we can see how the $\mathrm {GLC}$ (Fig. 10a) enters the cellular domain, diffuses and takes part in $\mathrm {HXK}$. On the other hand, the region where the reaction $\mathrm {act}$ (consuming $\mathrm {ATP}$ and producing $\mathrm {ADP}$) is clearly visible (Fig. 10b, c). Observing the plot of $\mathrm {GLY}$ (Fig. 10d) we can see the region where $\mathrm {GLY}$ is produced by $\mathrm {HXK}$ and at the same time consumed by $\mathrm {PYRK}$. Interestingly, we can notice that at the initial time, the $\mathrm {PYR}$ produced by $\mathrm {PYRK}$ contributes to generating $\mathrm {LAC}$, however $\mathrm {PYR}$ has not reached the $\mathrm {Mito}$ reaction site, so that there is no production of $\mathrm {ATP}$ from this reaction (Fig. 10e, f).

Fig. 10

Solutions of the Model obtained in a heart-shaped domain at the Initial time ($t=0$). $\mathrm {GLC}$ is entering the domain from the infimum of the domain. $\mathrm {ATP}$ is used in the region where $\mathrm {act}$ and $\mathrm {HXK}$ are located, vice versa $\mathrm {ADP}$ is produced. $\mathrm {GLY}$ is produced from $\mathrm {HXK}$. We can see $\mathrm {PYR}$ and $\mathrm {LAC}$ produced by $\mathrm {PYRK}$ and $\mathrm {LDH}$, respectively

Fig. 11

Solution of the Model obtained in a heart-shaped domain at the final time (T = 10). From top to bottom from left to right: $[\mathrm {GLC}],[\mathrm {ATP}],[\mathrm {ADP}],[\mathrm {GLY}],[\mathrm {PYR}],[\mathrm {LAC}]$. $\mathrm {GLC}$ is still consumed in the infimum where $\mathrm {HXK}$ is. $\mathrm {ADP}$, $\mathrm {PYR}$ and $\mathrm {LAC}$ have reached a state where the amount of concentration is uniformly distributed into the domain. $\mathrm {ATP}$ is consumed to produce $\mathrm {ADP}$ by $\mathrm {act}$ and produced by the $\mathrm {Mito}$ reaction. $\mathrm {GLY}$ is consumed by $\mathrm {PYRK}$ and produced by $\mathrm {HXK}$

At the final time ($T=10$), the $\mathrm {GLC}$ has diffused into the domain, and is still used to generate $\mathrm {HXK}$ (Fig. 11a). Concentrations $\mathrm {ADP}$, $\mathrm {PYR}$ and $\mathrm {LAC}$ have reached a stable state (Fig. 11c, e, f), while $\mathrm {ATP}$ is consumed by the $\mathrm {act}$ reaction and is produced by the $\mathrm {Mito}$ reaction (Fig. 11b). $\mathrm {GLY}$ is still consumed by $\mathrm {PYRK}$ and produced by $\mathrm {HXK}$ (Fig. 11d).

Briefly, we observe that CutFEM solves the test case inside the heart shaped domain, the singularities are well described thanks to the choice of a fine mesh. The computational time of CutFEM in this test case is smaller than in the previous experiment since the number of DOFS is 4 times less then the previous test case.

Experiment three: 3 dimensional complex domain

For the last experiment, we work with a 3 dimensional domain. We consider a union of six spheres that generate a pop-corn shape. To create the shape we use the union of the following level set functions

$$\begin{aligned} \phi _1(x,y,z)&= (x-4)^2 + (y-3)^2 + z^2 - 25\\ \phi _2(x,y,z)&= (x-12)^2 + (y-3)^2 + z^2 - 25\\ \phi _3(x,y,z)&= (x-8)^2 + (y-7)^2 + z^2 - 25\\ \phi _4(x,y,z)&= (x-8)^2 + (y+1)^2 + z^2 - 25\\ \phi _5(x,y,z)&= (x-8)^2 + (y-3)^2 + (z+4)^2 - 25\\ \phi _6(x,y,z)&= (x-8)^2 + (y-3)^2 + (z-4)^2 - 25 \end{aligned}$$

We sketch in Fig. 12a the problem setting projected onto the xy plane. The parameters chosen for the reaction locations are: $\mathcal {G}_{\mathrm {HXK}}(x_0=0.5, y0=2.0, z_0=0.0, \sigma = 0.5)$, $\mathcal {G}_{\mathrm {PYRK}}(x_0=1.1, y_0=1.2, z_0=0.0, \sigma =0.5)$, $\mathcal {G}_{\mathrm {LDH}}(x_0=4.0, y_0=5.0, z_0=0.0,\sigma = 0.5)$, $\mathcal {G}_{\mathrm {Mito}}(x_0=4.0,y_0=7.5, z_0=0.0, \sigma =0.5)$ and $\mathcal {G}_{\mathrm {act}}(x_0=6.0,y_0=8.5, z_0=0.0, \sigma =0.5)$. The source term is located in a ball of radius 0.3 centered on the origin, and we increase the source influx to 1000 until time 1. The diffusive constants, the reaction rates and the initial conditions are set as in “Asymptotic solution ODEs” section.

Fig. 12

a Setting of chemical reactions Eqs. (1)–(5) and the source term of $\mathrm {GLC}$ for the 3D domain, projected into the x–y axes. b The plot shows the average concentrations of each species computed in the three-dimensional example (blue line) compared with the one obtained in Fig. 6b with CutFEM in the circle

Based on the converged mesh sizes obtained for the two-dimensional examples, we set the mesh size to $h=0.3772$. Note that the choice of the Gaussian parameters $\sigma$ is designed to make the reactions significant since the mesh is coarser than the meshes used in the previous examples. Also, we increase the number of time steps to 100, while the final time is 10. The CutFEM penalty parameter is 0.1.

In Fig. 12b, we plot the average concentration of each species using Eq. (23). We examine the trend of the average concentrations plotting them with the ones computed in “Example one: two-dimensional circular domain” section. Even though the system solved one in 2D and one in 3D gives different results, we expect their average concentrations to behave similarly. Clearly, we have a higher influx of $\mathrm {GLC}$ entering the three-dimensional domain, as well as higher initial concentrations of $\mathrm {ATP}$ and $\mathrm {ADP}$. As expected, the average concentrations are higher in the three-dimensional domain but they all follow the same trend as the two-dimensional cases.

We show the results obtained in the three-dimensional domain at time 5 in Fig. 13. The results are in accordance with the ones showed in Section . We notice that the species are diffusing inside the domain and the reactions are happening in the location we indicated with the Gaussian functions.

Fig. 13

Solution of the reaction–diffusion system at time $t=5$ solved using CutFEM inside a 3 dimensional pop corn shape. From top to bottom from left to right $[\mathrm {GLC}]$, $[\mathrm {ATP}]$, $[\mathrm {ADP}]$, $[\mathrm {GLY}]$, $[\mathrm {PYR}]$, $[\mathrm {LAC}]$. The $\mathrm {GLC}$ is diffusing into the domain starting from the influx location. $\mathrm {ATP}$ and $\mathrm {ADP}$ are consumed and produced by the reaction $\mathrm {act}$. $\mathrm {GLY}$ is produced by $\mathrm {HXK}$ and $\mathrm {PYR}$ produced by $\mathrm {PYRK}$

In conclusion, in our three-dimensional experiment we can see how the chain of chemical reactions describe the energy metabolism of a cell in a complex structure. CutFEM is able to work with complex three-dimensional geometries, although the computational time is large (is 3 days in this case).

Conclusion and discussion

Studying brain energy metabolism can be helpful to further our understanding of several neurodegenerative diseases such as Alzheimer’s and Parkinson’s but a mechanistic perspective is still lacking. To investigate energy metabolism in biological relevant morphologies, we presented a simplified model of metabolism in a single cell using reaction–diffusion equations and demonstrated how to discretize the model using the FEM and CutFEM numerical methods. CutFEM  has the advantage of enabling the use of a single level set function per cell, independent of the complexity of the cell geometry. Implementing using FEniCS and libCutFEM ensures the usability by non-experts. In particular, modifying one or several of the reaction diffusion equations can be made by altering only a few lines of code and any linearisation and parallelisation is done automatically by the FEniCS framework. One of the difficulties of non-conforming methods for implicit domain definition [40, 61, 62] is due to the ill-conditioning of the system matrices when the interfaces or the boundary of the domain passes close to a node, leading to very large or very small diagonal terms. The consistent stabilisation term used by CutFEM prevents this issue and leads to a stable and convergent scheme. Other approaches are provided in [42,43,44].

Our results indicate that CutFEM is a valuable approach to deal with biological problems with arbitrarily complex cell morphologies. The appeal of level set descriptions coupled to enriched finite element approaches such as CutFEM lies in the fact that the motion of geometries over time only requires updating the level set function, without modifying the mesh. This will be instrumental to consider complex and time dependent shapes of cells such as astrocytes.

Our next steps will be to initiate the geometry of the cells through microscopy images 1 and to investigate the effect of perturbed astrocytic metabolism on neuronal support, which could lead to a better understanding of mechanisms in neurodegeneration.

Appendix A: Asymptotic solution

The system of ordinary differential equations is obtained applying the law of mass action to Eqs. (1)–(5). To avoid misinterpretation, we denote the concentration $A=[\mathrm {GLC}]$, $B =[\mathrm {ATP}]$, $C=[\mathrm {ADP}]$, $D=[\mathrm {GLY}]$, $E = [\mathrm {PYR}]$, $F=[\mathrm {LAC}]$ and we obtain the system

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{d A}{d t}=- K_{HXK}AB^2\\ \frac{d B}{d t}=-2K_{HXK} AB^2 +2K_{PYRK} C^2D + 28 K_{Mito}E- K_{act}B \\ \frac{d C}{d t}=2K_{HXK}AB^2 -2K_{PYRK}C^2D+K_{act}B \\ \frac{d D}{d t}=2K_{HXK}AB^2 -K_{PYRK}C^2D\\ \frac{d E}{d t}=K_{PYRK}C^2D -K_{LDH} E- K_{Mito}E\\ \frac{d F}{d t}= K_{LDH} E \end{array}\right. } \end{aligned}$$

(24)

The initial condition for Eq. (24) are chosen using the solution of the PDEs system Eq. (6). We compute the integral over the domain $\Omega$ using Eq. (23) of the solutions of the PDEs when the source term ends, making stable the total amount of concentrations inside the domain

$$\begin{aligned} {\left\{ \begin{array}{ll} A(t=0) =\int _\Omega [GLC](t=1)\text {d}x := A_0\\ B(t=0) =\int _\Omega [ATP](t=1)\text {d}x :=B_0\\ C(t=0) =\int _\Omega [ADP](t=1)\text {d}x := C_0\\ D(t=0) =\int _\Omega [GLY](t=1)\text {d}x :=D_0\\ E(t=0) =\int _\Omega [PYR](t=1)\text {d}x :=E_0\\ F(t=0) =\int _\Omega [LAC](t=1)\text {d}x := F_0 \end{array}\right. } \end{aligned}$$

(25)

Considering abundance of $\mathrm {ATP}$ and $\mathrm {ADP}$ inside the domain, this lead to only a possible system of steady state solution, that is when all the $\mathrm {GLC}$ is consumed by reaction $\mathrm {HXK}$ Eq. (1), $\mathrm {GLY}$ is consumed by reaction $\mathrm {PYRK}$ Eq. (2) and $\mathrm {ATP}$ is fully transformed into $\mathrm {ADP}$ from Eq. (5):

$$\begin{aligned} {\left\{ \begin{array}{ll} A_\infty = 0\\ B_\infty = 0\\ C_\infty = C_0 +B_0 + 28 (E_0+D_0+2A_0)\alpha \\ D_\infty = 0\\ E_\infty = 0 \\ F_\infty = F_0 + (E_0 + 2A_0 + D_0) (1-\alpha )\end{array}\right. } \end{aligned}$$

(26)

In the solutions we have introduced a parameter $\alpha \in [0,1]$ which take into account the fact that the concentration of $\mathrm {PYR}$ is transformed not in equal part to $\mathrm {LAC}$ and $\mathrm {ATP}$ from the reaction $\mathrm {LDH}$ Eq. (3) and reaction $\mathrm {Mito}$ Eq. (4).