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

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 develop- 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] ments 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.
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 reactiondiffusion 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.  [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 quasiindependently of the mesh (see Fig. 2).
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 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 (GLC) into two molecules of pyruvate (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 (LAC) by the enzyme lactate dehydrogenase (LDH) simplified with the chemical reaction (3) or enter mitochondria and used to produce ATP through the TCA cycle shown in reaction (4).
Last, we take into account the activity of the cell which uses 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. The main pathways we consider in our Simplified Model of Energy Metabolism of a cell. The glucose (GLC) enters inside the cell into the cytosol and takes part into the glycolysis composed by the two reactions HXK and PYRK. The products of the glycolysis, PYR and ATP, are then used into LDH reaction, that generates LAC, into act reaction, that describes the cellular activity producing ADP, and inside the mitochondria where the reaction denoted Mito happens generating the major source of ATP for the cell

Strong formulation of governing equations
Let be an open and bounded subset of R d (d = 2 or 3) with Lipschitz-continous boundary and we denote the concentration of the chemical species using the bracket notation, [ · ], as a function 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 [ · ] > 0 the diffusive constants. The reactions obey to the law of mass action [55] where K HXK , K PYRK , K LDH , K Mito and K act are the rate constants and we introduce Gaussian functions, indicated with G HXK (x 0 , σ ), G PYRK (x 0 , σ ), G LDH (x 0 , σ ), G Mito (x 0 , σ ) and G act (x 0 , σ ), to locate the region where the reactions are happening where x 0,· ∈ · and σ ∈ R. Eventually, to represent the entrance of the glucose inside the cytosol we define a source term function f : The strong form of the reaction diffusion system is then, finding the concentrations The system is completed with homogeneous Neumann boundary conditions on ∂ and initial conditions at time t = 0 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 [·] 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 ∈ H 1 ( ) and we denote with The weak form of system (6) can be found multiplying each equation by a test function Since D [·] is constant, we apply integration by parts to the second-order spatial derivative Since we have specified pure Neumann boundary conditions for each concentration species, the boundary terms vanish. Leading to We can, then, substitute Eq. (10) into Eq. (8) and use the following compact notation 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 t the size of the timestep. We denote with subscript n a concentration of (possibly multiple) chemical species at time t n , where 0 ≤ n < N with N = T / t is an integer that counts the time steps. The backward Euler method then approximates the continuous time derivative as We discretize the spatial domain using a triangulation T h that conforms (matches) the exact boundary ∂ . On this mesh we define the space of piece wise Lagrange finite elements of degree one as V and with V the product The space V can then be used to discretize the weak form equations (12). We The initial conditions are denoted [U ] 0 . We then solve a sequence of problems: where 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 using a mesh that conforms to the boundary, in CutFEM, the problem domain is described by a level-set function. The level set function is a scalar function on R d , such that φ(x) < 0 for x ∈ , φ(x) > 0 for x / ∈ and φ(x) = 0 for x ∈ ∂ . We then cover the domain by a regular background mesh ( ⊆ ) of simple shape, e.g. a box containing meshed with a uniform triangulation. Let K denote a triangle/tetrahedron in this triangulation. Now, letT h be the fictitious domain mesh composed by all elements K ∈ such that K ∩ = 0 ( ⊆T h ⊆ ). Furthermore, the union of all elements inT h is called the fictitious domain˜ . We denote with the set of elements intersected by the interface, and we define the set of so-called ghost penalty facets [46] (see Fig. 2) The stabilisation term introduced in the next section will be applied to this subset of facets. We consider the space We then solve a sequence of problems: inspired by [47]. Here, γ is a positive penalty parameter. These stabilization terms extend the solution from the physical domain onto the fictitious domain˜ , it is consistent with the continuous system. They prevent ill-conditioning of the system matrices in case only small parts of are contained in an element near the boundary ∂ . 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 Cut-FEM 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 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 spaceV 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 and in parts of the cut cells inside . 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 F G (Fig. 4). 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 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  (3) and (4), in order to obtain that PYR is used equally in the two reactions, such that we can set the parameter α of the ODE asymptotic solutions of Eq. (26) equal to 0.5.
The influx of 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 · of the chemical reactions are set equal to 10.0, in order to avoid that one reaction dominates over, and the diffusive parameters D [ · ] 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 [ATP] and [ADP] that are equal to 1. Final time is 1000, as well as the number of time step. The CutFEM penalty parameter γ 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 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 to obtain the average chemical concentration of each species at each time step. We use the following notatioñ 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 GLC, ATP, GLY and PYR go to zero whilst the products of the system are LAC and 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, , as a circle defined by the level set function φ(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 φ(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 Mito reaction that is located at G Mito (x 0 = 4.0, y 0 = 7.5, σ = 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 [ · ] = 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 σ .
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 GLC and 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 ATP is consumed in two regions that locate respectively HXK and act, while ADP is produced. The GLC entered the cell is diffusing in the domain and taking part to the HXK reaction. In comparison, the results of FEM and CutFEM are almost identical.
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.
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 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 φ(x, The description of the problem setting is shown in Fig. 9. We consider a source term for GLC active within a disk of radius 0.3, centered at the cusp and take the GLC influx to be 100. 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. The source term of GLC is no longer active, and the quantity of GLC is spreading through the domain and participating into the reaction HXK. ATP and ADP show that in the region where HXK, PYRK and act are located, when one is consumed the other is produced. We can notice production of GLY, PYR and LAC. Comparing the results in the three domain, they are extremely close to each other and visually not distinguishable  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 GLC (Fig. 10a) enters the cellular domain, diffuses and takes part in HXK. On the other hand, the region where the reaction act (consuming ATP and producing ADP) is clearly visible (Fig. 10b, c). Observing the plot of GLY (Fig. 10d) we can see the region where GLY is produced by HXK and at the same time consumed by PYRK. Interestingly, we can notice that at the initial time, the PYR produced by PYRK contributes to generating LAC, Fig. 9 Location of the chemical reactions Eqs. (1)-(5) and the source term of GLC in the heart shaped domain however PYR has not reached the Mito reaction site, so that there is no production of ATP from this reaction (Fig. 10e, f).
At the final time (T = 10), the GLC has diffused into the domain, and is still used to generate HXK (Fig. 11a). Concentrations ADP, PYR and LAC have reached a stable state (Fig. 11c, e, f), while ATP is consumed by the act reaction and is produced by the Mito reaction (Fig. 11b). GLY is still consumed by PYRK and produced by 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 10 Solutions of the Model obtained in a heart-shaped domain at the Initial time (t = 0). GLC is entering the domain from the infimum of the domain. ATP is used in the region where act and HXK are located, vice versa ADP is produced. GLY is produced from HXK. We can see PYR and LAC produced by PYRK and LDH, respectively  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.
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 σ 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 GLC entering the three-dimensional domain, as well as higher initial concentrations of ATP and 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.
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 reactiondiffusion 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. . The GLC is diffusing into the domain starting from the influx location. ATP and ADP are consumed and produced by the reaction act. GLY is produced by HXK and PYR produced by PYRK 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.

Abbreviations
Considering abundance of ATP and ADP inside the domain, this lead to only a possible system of steady state solution, that is when all the GLC is consumed by reaction HXK Eq.
(1), GLY is consumed by reaction PYRK Eq. (2) and ATP is fully transformed into ADP from Eq. (5): In the solutions we have introduced a parameter α ∈ [0, 1] which take into account the fact that the concentration of PYR is transformed not in equal part to LAC and ATP from the reaction LDH Eq. (3) and reaction Mito Eq. (4).