Claire Guerrier
email: claire.guerrier@univ-cotedazur.fr

Stella Krell

Paul Paragot

Nonlinear coupling of the Poisson-Nernst Planck system of equations using the Discrete Duality Finite Volumes method, application to ionic and voltage dynamics in neuronal compartments

Keywords: Poisson-Nernst Planck, Discrete Duality Finite Volumes, Dendritic integration, Neuronal plasticity

come

Nonlinear coupling of the Poisson-Nernst Planck system of equations using the Discrete Duality Finite Volumes method, application to ionic and voltage dynamics in neuronal compartments

Introduction

The Poisson-Nernst Planck (PNP) system of equations is a well-established framework for describing the movement of ions and the variations of the electric field in an electrolyte [START_REF] Kirby | Micro-and nanoscale fluid mechanics: transport in microfluidic devices[END_REF]. It is employed in various fields such as semiconductor devices [START_REF] Jüngel | Quasi-hydrodynamic semiconductor equations[END_REF], electrophoresis models [START_REF] As Khair | Ion steric effects on electrophoresis of a colloidal particle[END_REF], or biological systems [Bazant et al., 2004, Eisenberg and[START_REF] Eisenberg | Poisson-Nernst-Planck systems for ion channels with permanent charges[END_REF]. Solving the PNP system analytically and numerically is challenging due to the nonlinear coupling between the electrostatic potential and the concentrations of ionic species. Also, the system can present, depending on the domain and boundary conditions, a stiff boundary layer. A standard example is a binary electrolyte between parallel-plate electrodes, a system showing a boundary layer near the electrodes, characterized by significant ionic concentration gradients [START_REF] Mz Bazant | Diffuse-charge dynamics in electrochemical systems[END_REF].

In neuroscience, it is rather recent that models using the PNP system of equations are developed. They follow the latest innovations in experimental techniques now reaching the nano-scale, which create the need for modeling at the same precision [Savtchenko et al., 2017, Holcman and[START_REF] Holcman | The new nanophysiology: regulation of ionic flow in neuronal subcompartments[END_REF]. Indeed, the classical model describing voltage propagation in neurons is the Cable theory, a one-dimensional macroscopic model that makes an analogy between the neuron and a wire, assuming voltage propagates as in an RC electrical circuit [START_REF] Koch | Cable theory in neurons with active, linearized membranes[END_REF]. In particular, this model neglects the variations of ionic concentrations within the neuronal cytoplasm. Such an assumption is valid while considering large axons, but is not relevant in smaller neuronal structures [START_REF] Lp Savtchenko | Electrodiffusion phenomena in neuroscience: a neglected companion[END_REF], where we expect the influence of the ionic dynamics on the electric field to be important. In contrast to this one-dimensional approach, the PNP system of equations describes the dynamics of ionic charges due to both diffusion and electric field and is highly suitable to capture such dynamics in small neuronal compartments. Computational models of ionic electrodiffusion in neurons usually combine the PNP equations and Hodgkin-Huxley formalism that represent voltage-gated ionic channels [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. For example, [START_REF] Cl Lopreore | Computational modeling of three-dimensional electrodiffusion in biological systems: application to the node of Ranvier[END_REF] introduced such a model, to simulate the electrodiffusion of ions within a node of Ranvier during an action potential. The numerical study focused on an axon model with a diameter of 15 µm, utilizing classical finite volume techniques with Delaunay-Voronoi dual meshes. The authors compared their simulations with a one-dimensional model based on the Cable equations coupled to Hodgkin-Huxley dynamics. Their simulation results showed remarkable similarities between the two models, emphasizing that action potential dynamics in axons are well described by a one-dimensional Cable model. [START_REF] Pods | Electrodiffusion models of neurons and extracellular space using the Poisson-Nernst Planck equations-numerical simulation of the intra-and extracellular potential for an axon model[END_REF] extended this exploration by focusing on studying the boundary layer (BL) formed near the membrane of neuronal axons. This two-dimensional investigation explicitly addressed BL phenomena using finite element methods, with the Hodgkin-Huxley model applied as a Neumann boundary condition on the membrane. They observed that the dynamics are very similar to the Cable model, except for some deviations close to the membrane, most probably caused by the influence of the intracellular potential on the extracellular space. These works focused on action potential dynamics in axons but did not address the case of low amplitude ionic influx in smaller neuronal compartments such as dendrites or dendritic spines. As experimental investigations of voltage and ionic dynamics in these small neuronal compartments remain challenging, mainly because of their microscopic scale, it is crucial to develop models dedicated to their specific dynamics, which is the aim of the present paper.

Many numerical methods have been developed to simulate the PNP equations in the past decades, focusing on various properties and applications: [START_REF] Bessemoulin-Chatard | Numerical analysis of ddfv schemes for semiconductors energy-transport models[END_REF][START_REF] Cao | An adaptive conservative finite volume method for Poisson-Nernst Planck equations on a moving mesh[END_REF][START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF][START_REF] Mathur | A multigrid method for the Poisson-Nernst Planck equations[END_REF][START_REF] Mirzadeh | A conservative discretization of the Poisson-Nernst Planck equations on adaptive cartesian grids[END_REF][START_REF] Moatti | A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes[END_REF][START_REF] Shen | A decoupling two-grid method for the time-dependent Poisson-Nernst-Planck equations[END_REF][START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF][START_REF] Su | A positivity-preserving and free energy dissipative hybrid scheme for the poisson-nernst-planck equations on polygonal and polyhedral meshes[END_REF][START_REF] Zheng | Second-order Poisson-Nernst Planck solver for ion transport[END_REF]. For our specific problem of voltage and concentration dynamics in neurons, we develop in this study a finite volume approach to solve the PNP system of equations, using the Discrete-Duality Finite Volume (DDFV) framework, and preserving the positivity of ionic concentrations. The DDFV method, introduced by [START_REF] Hermeline | A finite volume method for the approximation of diffusion operators on distorted meshes[END_REF] and [START_REF] Domelevo | A finite volume method for the laplace equation on almost arbitrary two-dimensional grids[END_REF], enables the study of the Laplace equation on various 2D meshes, accommodating nonconformal and deformed meshes. Indeed, DDFV schemes double the discrete unknowns, setting unknowns at both the vertices and the centers of the initial mesh. This set of unknowns allows for the definition of a full discrete gradient on a new mesh called the diamond mesh. Hence, the method can be used on general meshes that do not necessarily satisfy classical orthogonality conditions. The DDFV methods construct two discrete operators, the gradient and the divergence, that are in duality in a discrete sense as in the continuous setting [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lionstype elliptic problems on general 2d meshes[END_REF]. Several studies have demonstrated the numerical advantages of DDFV methods for PDEs, and these methods have been developed for a large number of models. The DDFV method is particularly attractive for numerically solving the PNP equations on a complex geometry, because of its local conservation property, and robustness with respect to mesh distortion. A telling example concerns a DDFV scheme discretizing the Peaceman model -that describes miscible displacement in porous media -on which a convergence result is proved [START_REF] Chainais-Hillairet | Convergence analysis of a ddfv scheme for a system describing miscible fluid flows in porous media[END_REF]. Numerical analysis of systems similar to PNP in the field of semiconductors was carried out by the authors [START_REF] Bessemoulin-Chatard | Numerical analysis of ddfv schemes for semiconductors energy-transport models[END_REF]. Then, a positivity-preserving scheme was introduced in [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing discrete duality finite volume scheme for convection diffusion equations[END_REF], to investigate the existence and large-time behavior of solutions in Fokker-Planck equations. Positivity was ensured by using a non-linear reformulation of the equation.

In this paper, we numerically solve the PNP system of equations using a DDFV approach, with a scheme preserving the positivity of ionic concentrations, and a high level of mesh refinement at the boundary layer when necessary. We challenge the accuracy of our DDFV scheme through several test cases and obtain numerically a second-order accuracy in space. We investigate the convergence and robustness of our method by comparing our results with a test case proposed in [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF] with a 2D finite difference scheme. We then compare our results with the 1D finite difference scheme of [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF] in the case of the presence of a boundary layer. Finally, we study voltage and ionic concentration dynamics within dendrites and dendritic spines. We first consider the propagation of an influx of ions at a dendritic bifurcation. We determine the distance between the bifurcation and an ionic reservoir such that the latter has a small influence on the dynamics, which is around 20 microns for our geometry. We then show that two signals arriving simultaneously at the bifurcation sum linearly. In a second time, we consider dendritic spines and investigate the influence of spine geometry on voltage and ionic concentration dynamics. We observe that if a 1D approximation for both voltage and concentration is relevant in the neck, the dynamics in the head require at least a 2D approach. We finally investigate the influence of two neighboring spines on each other, when an influx of ions arrives in one spine. We observe that depending on the distance to an ionic reservoir, a dendritic spine can either act as an autonomous compartment or be subject to signal invasion from the neighboring spine.

The paper is organized as follows: Section 2 presents the formulation for the PNP system of equations in the case of two ion species. Section 3 provides a detailed description of the numerical scheme. Section 4 presents the numerical simulations and convergence results. Section 5 focuses on applications in neuroscience. Finally, the discussion and conclusion are presented in the last section.

Presentation of the Poisson-Nernst Planck system of equations

In this section, we present the PNP system of equations, which describes the movement of ions in an electrolyte. The Nernst Planck equation links the distribution of the ionic concentrations c i of species i to the electric potential V through the ion conservation equation, using a drift-diffusion flux J i . J i divides into a diffusion term that models the Brownian motion of ions according to Fick's law, as well as a convective term that describes the ion transport resulting from the electric field. This set of equations is then combined with the Poisson equations for voltage dynamics, creating the PNP system of equations. It is worth noting that the PNP system consists of a combination of parabolic and elliptic equations. We consider in this paper the case of two ionic species P and N , with respective valences +1 and -1, and concentration c P and c N -our method can be extended to any number of ions without difficulties.

Let Ω ⊂ R 2 be a connected open bounded domain and T f > 0 a finite time. The PNP system of equations on Ω × [0, T f ] writes as:

∂ t c P = -∇ • J c P + f c P , ∂ t c N = -∇ • J c N + f c N , -∇ • (εε 0 ∇V ) = F (c P -c N ) + f V , (1)
with:

J c P = -D P ∇c P - F D P RT θ c P ∇V, J c N = -D N ∇c N + F D N RT θ c N ∇V.
(2)

In system (1), the first (resp. second) equation represents the dynamics of species P (resp. N ), with drift-diffusion flux J c P (resp. J c N ) and source term f c P (resp. f c N ). The third equation corresponds to the Poisson equation, where ε 0 is the vacuum permittivity, ε the relative permittivity of the solution, and f V a source term. Eqs. (2) describe the drift-diffusion fluxes, with D P > 0 (resp. D N >0) the diffusion coefficient of species P (resp. species N ), R the gas constant, T θ the absolute temperature, and F the Faraday constant. In the rest of the paper, we will consider that the source terms

f c P , f c N , f V are in L 2 (Ω × T f ).
We define:

γ = ϵϵ 0 RT θ F 2 , β = F RT θ , (3)
such that the system rewrites:

∂ t c P = -∇ • (D P (∇c P + c P β∇V )) + f c P , ∂ t c N = -∇ • (D N (∇c N -c N β∇V )) + f c N , -∇ • (γβ∇V ) = c P -c N + f V . (4)
We use a positivity-preserving DDFV scheme as presented in [START_REF] Cancès | Numerical analysis of a nonlinear free-energy diminishing discrete duality finite volume scheme for convection diffusion equations[END_REF] to ensure the positivity of ionic concentrations. The method is based on the reformulation of the fluxes J c P and J c N :

J c P = -D P c P ∇ (log c P + βV ) , J c N = -D N c N ∇ (log c N -βV ) .
(5) Finally, the system (1), defined on Ω × [0, T f ], rewrites:

∂ t c P = ∇ • (D P c P ∇ (log c P + βV )) + f c P , ∂ t c N = ∇ • (D N c N ∇ (log c N -βV )) + f c N , -∇ • (γβ∇V ) = c P -c N + f V .
(6)

Initial conditions and boundary values

The initial conditions of the PNP system (6), defined on Ω, are:

c P (t = 0, x) = c 0 P (x), c N (t = 0, x) = c 0 N (x), V (t = 0, x) = V 0 (x) ∀x ∈ Ω, (7)
with

c 0 P , c 0 N , V 0 ∈ L 2 (Ω).
We decompose the domain boundaries as ∂Ω = Γ Dir ∪Γ N eu , where Γ Dir and Γ N eu represent respectively the subsets of the boundary corresponding to Dirichlet and Neumann conditions. The Dirichlet boundary conditions write as:

c P = c Dir P , c N = c Dir N , V = V Dir , (8) with c Dir P , c Dir N , V Dir ∈ L 2 (Γ Dir × [0, T f ]).
Neumann boundary conditions write as:

J c P • n = g, J c N • n = h, ∇V • n = 0, (9)
where n is the unit outward normal to Ω, and g, h

∈ L 2 (Γ N eu × [0, T f ]).
3 Presentation of the scheme

Meshes and notations

In this section, we present the DDFV scheme and introduce three distinct meshes: the primal mesh, the dual mesh, and the diamond mesh (see Figure 1 and 2). The primal mesh, denoted as M = M ∪ ∂M, is composed of the interior primal mesh M and the boundary primal mesh ∂M. The interior primal mesh M is a partition of Ω with polygonal cells. The boundary primal mesh ∂M is the set of boundary edges of M. Note that these cells are treated as degenerate cells.

For each element K ∈ M, we denote x K as the center of the cell K. The set of vertices of the primal mesh is denoted as X * . We differentiate between interior vertices and boundary vertices: for a given x K * ∈ X * such that x K * is not on the boundary ∂Ω, we define a dual cell K * by connecting the centers of all K that share x K * as a vertex. This set of polygons K * forms the interior dual mesh M * . For each x K * ∈ X * ∩ ∂Ω, i.e. x K * is on the boundary, we construct a dual cell K * by joining x K * and the x K where K ∈ M ∪ ∂M, that share x K * as a vertex. This set of polygons defines the boundary dual mesh ∂M * . We denote the dual mesh as M * = M * ∪∂M * . Fig. 1 illustrates two examples of a primal and a dual mesh. The measure of cell K (resp. K * ) is denoted as

|K| (resp. |K * |).
For two neighboring primal cells, K and L, we consider the intersection of their boundaries, denoted as ∂K ∩ ∂L, which corresponds to a line segment and is referred to as an edge of the mesh M. We denote this edge as σ = K|L. The set of edges that lie within the interior of

• • • • • • • • • • • • • • • • • • • x K x L Primal mesh M = M ∪ ∂M ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ ⊠ • • • • • • • x K x L • • • • • • • • • • • • x K * x L * Dual mesh M * = M * ∪ ∂M * Figure 1: Left:
Example of a primal mesh (grey). Black dots are the centers of each primal cell K, denoted by x K . One cell is colored in light grey. Right: the corresponding dual mesh (blue). Blue dots are the centers of each dual cell K * (resp. L * ) which are denoted by x K * (resp. x L * ). An example of a dual cell is colored in light blue.

the domain is denoted as E int . The set of primal edges on the boundary is denoted E ext . The total set of primal edges is E = E int ∪ E ext . Similarly, the set of edges in the interior of the dual mesh is denoted as E * int , the set of edges on the boundary of the dual mesh is E * ext , and we define E * = E * int ∪ E * ext , the total set of dual edges. For each pair (σ, σ * ) ∈ E × E * , where σ = (x K * , x L * ) and σ * = (x K , x L ), we define a quadrilateral diamond cell, denoted as D, with σ and σ * as its diagonals. In the case where σ is a boundary edge, the diamond D degenerates into a triangle. The entire set of diamond cells forms the diamond mesh, denoted as D, which can be expressed as D = D ext ∪ D int . Here, D ext represents the set of boundary diamonds, where σ ⊂ ∂Ω, and D int represents the set of interior diamonds, where σ ̸ ⊂ ∂Ω. In summary, the DDFV mesh T consists in the union of the primal mesh M and the dual mesh M * , along with the diamond mesh D (Fig. 1). For a diamond D with vertices (x K * , x L , x L * , x K ), we define x D as its center, |D| as its measure, |σ| as the length of the primal edge σ, |σ * | as the length of the dual edge σ * , and α D as the angle between the vectors (x L , x K ) and (x L * , x K * ). Additionally, we define two mutually orthogonal basis vectors, (τ K * ,L * , n σK ), and (n σ * K * , τ K,L ), where n σK is the unit vector normal to σ oriented from K to L, n σ * K * is the unit vector normal to σ * oriented from K * to L * , τ K * ,L * is the unit tangent vector to σ oriented from x K * to x L * , and τ K,L is the unit tangent vector to σ * oriented from x K to x L (Fig. 2). It is worth noting that the area of a diamond cell D is given by |D| = 1 2 |σ||σ * | sin(α D ). Furthermore, we define two subsets for the diamond cells as follows:

• • • • • • • • • • • • • • • • • • • x K x L x K * x L * • • • • x L * x K * x L x K τ K * ,L * σ n σK σ * τ K,L n σ * K *
D K = {D ∈ D, σ ∈ E K }, D K * = {D ∈ D, σ * ∈ E K * }.
where:

-E K is the set of edges of K ∈ M ∪ ∂M. Note that for all K ∈ ∂M, E K is reduced to one edge, i.e. σ = K. -E K * is the set of edges of K * ∈ M * ∪ ∂M * .
We consider mixed boundary conditions, and thus define four subsets of the boundary mesh. The boundary primal and dual meshes for the Dirichlet boundary condition are:

∂M Dir = {K ∈ ∂M : x K ∈ Γ Dir }, ∂M * Dir = {K * ∈ ∂M * : x K * ∈ Γ Dir }.
The boundary dual mesh for the Neumann boundary condition is:

∂M * N eu = {K * ∈ ∂M * : x K * ∈ Γ N eu \Γ Dir }.
Finally, the boundary diamond mesh for the Neumann boundary condition is:

D ext,N eu = {D ∈ D ext : σ ∈ D ∩ Γ N eu }.
Note that per definition, there is a one-to-one relationship between the edges of primal cells and diamonds, and between the edges of dual cells and diamonds. We finally denote Z * S the set of dual cells on a segment S defined as:

Z * S = {K * ∈ M * , such that K * ∩ S ̸ = ∅} (10)
and K * X the unique dual cell containing a point X:

K * X = {K * ∈ M * , such that K * ∩ X ̸ = ∅} (11)

Discrete unknowns and operators

The DDFV method enables the construction of two-dimensional discrete gradient and divergence operators that exhibit duality in a discrete sense. A description of the duality framework can be found in [START_REF] Domelevo | A finite volume method for the laplace equation on almost arbitrary two-dimensional grids[END_REF] and [START_REF] Andreianov | Discrete duality finite volume schemes for Leray-Lionstype elliptic problems on general 2d meshes[END_REF]. In this subsection, we introduce the discrete unknowns, the discrete gradient and divergence operators, as well as the reconstruction operator, that maps the unknowns on the primal and dual meshes on the diamond mesh. We first describe the sets of discrete unknowns. We define R T the linear space of scalar fields that are constant over the cells of M and M * :

u T ∈ R T ⇐⇒ u T = (u K ) K∈M , (u K * ) K * ∈M * .
We set R 2 D the linear space of vector fields constant on the diamonds:

ξ D ∈ R 2 D ⇐⇒ ξ D = ξ D D∈D .
We also denote R D as the set of scalar fields that are constant over the diamonds. The discrete gradient operator is defined as a mapping from R T to R 2 D , denoted by ∇ D u T = (∇ D u T ) D∈D for all u T ∈ R T . For each diamond D ∈ D:

∇ D u T = 1 2|D| [|σ|(u L -u K )n σK + |σ * |(u L * -u K * )n σ * K * ]. (12)
The discrete divergence operator is a mapping from R 2 D to R T given by, for all ξ D ∈ R 2 D :

div T ξ D = div M ξ D , div ∂M ξ D , div M * ξ D , div ∂M * ξ D . ( 13
)
The divergence on the primal mesh is div M ξ D = (div K ξ D ) K∈M and div ∂M ξ D = 0, with:

div K ξ D = 1 |K| D∈D K |σ|ξ D • n σK , ∀K ∈ M.
The divergence on the dual mesh is as well div

M * ξ D = (div K * ξ D ) K * ∈M * and div ∂M * ξ D = (div K * ξ D ) K * ∈∂M * with: div K * ξ D = 1 |K * | D∈D K * |σ * |ξ D • n σ * K * , ∀K * ∈ M * , div K * ξ D = 1 |K * |   D∈D K * |σ * |ξ D • n σ * K * + D∈D K * ∩Dext |σ| 2 ξ D • n σK   , ∀K * ∈ ∂M * .
We introduce r D , a reconstruction operator on diamonds that maps R T to R D . For any

u T ∈ R T , r D [u T ] is defined as (r D (u T )) D∈D
, where D ∈ D and its vertices are denoted as

(x K , x K * , x L , x L * ): r D (u T ) = 1 4 u K + u L + u K * + u L * . ( 14
)

Discrete notations

Let N be a positive integer. We denote dt = T f /N and t n = n × dt for n ∈ 0, ..., N . We define the discrete space R T ,dt := R T N +1 and its associated discrete vector u T ,dt ∈ R T ,dt . Furthermore, for n ∈ N and a given function u, we define the discrete projection P T ,n m u = (P M,n m u, P M * ,n m u, P ∂Ω,n m u) as follows:

P M,n m u = 1 |K| K u(t n , x)dx K∈M , P M * ,n m u = 1 |K * | K * u(t n , x)dx K * ∈M * , P ∂Ω,n m u = 1 |K| K u(t n , x)dx K∈∂M , 1 |K * | K * u(t n , x)dx K * ∈∂M * .
Also, we give Dirichlet discrete projections:

P ∂M Dir ,n m u = 1 |K| K u(t n , x)dx K∈∂M Dir , P ∂M * Dir ,n m u = 1 |K * | K * u(t n , x)dx K * ∈∂M * Dir .

DDFV scheme for the PNP system

In this section, we present the DDFV scheme. We first discretize the source terms:

f T ,dt c P = P T ,dt m f c P , f T ,dt c N = P T ,dt m f c N , f T ,dt V = P T ,dt m f V ,
and the initial conditions:

c T ,0 P = P T ,0 m c 0 P , c T ,0 N = P T ,0 m c 0 N , V T ,0 = P T ,0 m V 0 .
Next, we define the discretization for the Dirichlet boundary condition functions: The same way, we define the discretization of the Neumann boundary condition functions g:

g Dext,dt N eu = (g σ,n ) σ∈D ext,N eu n∈{0,...,N } with g σ,n = σ g(t n , x)dx.
We now present the numerical scheme for the PNP system (6)-( 9). In order to ensure stability, we adopt an implicit Euler scheme for time discretization. At each time step, we employ the Newton method to solve the nonlinear system of equations. By integrating equations ( 6)-( 9) over M and M * ∪ ∂M * N eu , we find the solutions (c T ,dt P , c T ,dt N , V T ,dt ) ∈ R T ,dt that satisfy the following nonlinear problem for all n ≥ 0:

c K,n+1 P -c K,n P dt + div K J D,n+1 c P = f K,n+1 c P , ∀K ∈ M, (15)
c K * ,n+1 P -c K * ,n P dt + div K * J D,n+1 c P = f K * ,n+1 c P , ∀K * ∈ M * ∪ ∂M * N eu , (16)
c K,n+1 N -c K,n N dt + div K J D,n+1 c N = f K,n+1 c N , ∀K ∈ M, (17)
c K * ,n+1 N -c K * ,n N dt + div K * J D,n+1 c N = f K * ,n+1 c N , ∀K * ∈ M * ∪ ∂M * N eu , (18)
-div K (γβ∇ D V K,n+1 ) + c K,n+1 N -c K,n+1 P = f K,n+1 V , ∀K ∈ M, (19)
-div K * (γβ∇ D V K * ,n+1 ) + c K * ,n+1 N -c K * ,n+1 P = f K * ,n+1 V , ∀K * ∈ M * ∪ ∂M * N eu , (20)
J D,n+1 c P = -D P r D (c T ,n+1 P )∇ D log c T ,n+1 P + βV T ,n+1 , (21)
J D,n+1 c N = -D N r D (c T ,n+1 N )∇ D log c T ,n+1 N -βV T ,n+1 . (22)
The discrete mixed boundary conditions can be expressed as:

V K,n+1 = V K,n+1 Dir , c K,n+1 P = c K,n+1 P,Dir , c K,n+1 N = c K,n+1 N,Dir , ∀K ∈ ∂M Dir (23) V K * ,n+1 = V K * ,n+1 Dir , c K * ,n+1 P = c K * ,n+1 P,Dir , c K * ,n+1 N = c K * ,n+1 N,Dir , ∀K * ∈ ∂M * Dir (24) |σ|∇ D V T ,n+1 • n σK = 0, |σ|J D,n+1 c P • n σK = g σ,n+1 , |σ|J D,n+1 c N • n σK = 0, ∀D ∈ D ext,N eu . ( 25
)
The existence of the discrete solution of ( 15)-( 25) is an oncoming work, that is not the focus of the current paper. It requires a paper in itself, dedicated to the study of existence, positivity of solution and convergence, that will be inspired by [START_REF] Moatti | A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes[END_REF].

Performance of the DDFV scheme

In this section, we present numerical simulations to test the performance of our DDFV scheme ( 15)-( 25). In subsection 4.1, we evaluate the accuracy of the DDFV scheme by considering two test cases with known exact solutions of the PNP system ( 6)-( 9). In the first test case, we compare our numerical results with an analytical solution of the PNP equations. In the second one, we compare our numerical results with existing analytical and numerical results coming from [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF]. In subsection 4.2, we show the numerical behavior of our DDFV results in the presence of boundary layer, using a test case inspired by the one-dimensional one in [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF]. The coefficients β, D P and D N are set to 1 in all the section.

Convergence results

To evaluate the accuracy of our DDFV scheme ( 15)-( 25), we compare it with exact solutions at various mesh resolutions, using the discrete norm

L ∞ ((0, T f ); L 2 (Ω)), ∀u T ,dt ∈ R T ,dt : ∥u T ,dt ∥ T ,∞ = max n∈{0,...,N }   1 2 K∈M |K||u K,n | 2 + 1 2 K * ∈M * |K * ||u K * ,n | 2   1 2 .
The computational domain Ω is set to ]0, 1[ 2 and the final time T f = 1. We consider only Dirichlet boundary conditions on ∂Ω, i.e., Γ N eu = ∅, and we set γ = 1. For the first test case, we consider c P,ex , c N,ex and V ex , the exact solution of ( 15)-( 25), defined as follows:

           c P,ex (t, x, y) = 7x + 5 + 3t 2 , c N,ex (t, x, y) = x + 1 + t 2 , V ex (t, x, y) = -x 3 -2y 2 -t 2 x 2 . ( 26
)
with boundary conditions, initial conditions, and source terms being:

           f c P (t, x, y) = 6t 4 + (46x + 22)t 2 + 6t + 63x 2 + 58x + 20, f c N (t, x, y) = -2t 4 + (-10x -6)t 2 + 2t -9x 2 -10x -4, f V (t, x, y) = 0.
The projection on R T ,dt of these exact solutions, {c T ,dt P,ex , c T ,dt N,ex , V T ,dt ex } are defined by: c T ,dt P,ex = P T ,dt m c P,ex , c T ,dt N,ex = P T ,dt m c N,ex , V T ,dt ex = P T ,dt m V ex . We simulate the system ( 15)-( 25) on several increasingly refined meshes generated by the quadrangle mesh family. One example mesh is depicted in Figure 3, with a mesh size h = 2.341×10 -1 . The respective mesh size h for each mesh ranges from 4.692×10 -1 to 6.253×10 -2 . The time step dt is initially set to 10 -2 for the mesh with the largest h = 4.692 × 10 -1 . For subsequent meshes, the time step is divided by 4. Recall that at each time step, we run a Newton algorithm to solve our non-linear problem. This Newton algorithm converges, with on average 2 to 1 iterations depending on the mesh refinement. The maximal number of iterations for all the time steps over all the meshes is 3. The errors e c P , e c N and e V , as well as the convergence order for the first test case (26) on the quadrangle mesh family.

V Order 4.692E-1 1.0E -2 2.593E-3 • 1.443E-3 • 2.070E-2 • 2.341E-1 2.5E -3 8.109E-
To illustrate the convergence of our method, we define the error according to the norm L ∞ ((0, T f ); L 2 (Ω)), between a numerical vector y and the corresponding exact solution y ex :

e y = ∥y T ,dt -y T ,dt ex ∥ T ,∞ . (27)
The errors e c P , e c N , and e V , for the different meshes, as well as the convergence order obtained numerically are given in Table 1. For both the concentrations c P and c N and the electrostatic potential V , we observe a convergence of order two.

We then investigate the convergence and robustness of our method by comparing our results with a test case proposed in [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF]. In this paper, the authors introduced a finite difference scheme for approximating the PNP solutions in 2D while preserving the positivity of the ionic concentrations. The computational domain Ω is set to ] -1, 1[ 2 , the final time T f = 0.1 and γ = 1. The boundary conditions are all Dirichlet (i.e., Γ N eu = ∅). The exact solutions of the second test case are defined as follows:

           c P,ex (t, x, y) = e -t cos(2πx) sin(2πy) + 2, c N,ex (t, x, y) = e -t sin(2πx) cos(2πy) + 2, V ex (t, x, y)
= e -t sin(2πx) sin(2πy).

(28)

The boundary and the initial conditions, as well as the source terms are:

                   f c P (t, x, y) = -16 cos(2πx)e -t e -t π 2 cos(2πy) 2 - 3 4 sin(2πx) + sin(2πy) - π 2 2 + 1 16 , f c N (t, x, y) = cos(2πy)(8π 2 -1)e -t sin(2πx) + 16 cos(2πx) 2 - 3 4
π 2 e -2t sin(2πy) cos(2πy), f V (t, x, y) = e -t 8π 2 sin(2πy) + cos(2πy) sin( 2πx)e -t cos(2πx) sin(2πy).

We run simulations using the same family of Cartesian meshes and the same time step dt = h 2 as in [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF]. At each time step, our Newton algorithm converges with a maximum of 3 iterations. The average number of iterations is between 1 and 2, depending on the mesh refinement.

h dt e c P Order e c N Order e V Order 1E-1 1.0E-2 4.086E-3 • 4.077E-3 • 5.453E-3 • 5E-2 2.5E-3 1.038E-3 1.98 1.036E-3 1.98 1.372E-3 1.99 2.5E-2 6.25E-4 2.606E-4 1.99
2.600E-4 1.99 3.436E-4 2.00 1.25E-2 1.56E-4 6.519E-5 2.00 6.504E-5 2.00 8.593E-5 2.00

Table 2: The errors e c P , e c N and e V , as well as the convergence order for the second test case (28), using Cartesian meshes.

In Table 2, we give the value of the errors e c P , e c N , e V (eq. ( 27)). We obtain again a convergence of order 2 for c P , c N and V . For each mesh size, the L ∞ ((0, T f ); L 2 (Ω)) errorsbetween the numerical solutions of ( 15)-( 25) and the exact solution -obtained from our DDFV scheme is nearly 10 times smaller than the one obtained in [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF]. As an example, we obtain a value for e c P of 6.519 × 10 -5 to compare with a value of 3.093 × 10 -4 in [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF] (h = 1.25 × 10 -2 ). Put differently, to obtain an error with the same order of magnitude, our DDFV scheme requires a mesh with half the number of cells than the finite difference scheme proposed in [START_REF] Liu | A positivity-preserving, energy stable and convergent numerical scheme for the Poisson-Nernst Planck system[END_REF], which significantly reduces computation costs.

Comparative analysis of the DDFV scheme at the boundary layer

We consider in this part a test case with a boundary layer initially presented in [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF]. This test case is defined on the one-dimensional domain ]0, 1[, with the following boundary conditions:

c P (t, 0) = 1 + t, c N (t, 0) = 1, V (t, 0) = 0, c P (t, 1) = 1, c N (t, 1) = 1 + t, V (t, 1) = 0.
This system has two boundary layer regions, ]0, 0.01[ and ]0.99, 1[. In [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF], the authors propose a finite difference scheme preserving the positivity of the concentration. They adopt a mesh size ∆x = 8 × 10 -4 in the boundary layer region, and a larger one ∆x = 3 × 10 -3 in the bulk region ]0.01, 0.99[.

We build our third test case by adapting this 1D test case to our 2D framework. Our computational domain Ω is ]0, 1[ 2 , and we adapt the Dirichlet boundary conditions so that they also depend on the space variable:

c Dir P (t, x, y) = 1 + (1 -x)t, c Dir N (t, x, y) = 1 + xt, V Dir (t, x, y) = 0. ( 29
)
We set the source terms of the system (6) to zero. The initial conditions are computed using eq. ( 29) at t = 0. The final time T f = 1, with a time step dt = 10 -2 and γ = 0.01. We use a fixed non-uniform mesh with a mesh size in the x-direction of ∆x = 10 -3 in the boundary layer region, namely ]0, 0.01[×]0, 1[ and ]0.99, 1[×]0, 1[, gradually increasing to ∆x = 10 -2 in the bulk region ]0.01, 0.99[×]0, 1[. For all cells, we set the mesh size in the y-direction equal to ∆y = 10 -2 . The mesh has 11800 quadrilateral cells (rectangles). At each time step, the Newton algorithm converges with a maximum of 2 iterations.

In Figure 4, we plot on the left (A-C-E) the dual values of the concentrations c P and c N and the electrostatic potential V at final time ) K * ∈Z * y=0.5 ), as well as the dual discrete electrostatic potential (V K * ,50 ) K * ∈Z * y=0.5 (resp. (V K * ,100 ) K * ∈Z * y=0.5 ) at t = 0.5 (resp. t = 1) with a linear trend of the potential. We recall that the definition of Z * y=0.5 is given by eq. ( 10). This figure can be compared to Fig. 2 in [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF].

T f = 1, that is (c K * ,100 P ) K * ∈M * , (c K * ,100 N ) K * ∈M * and (V K * ,100 ) K * ∈M * .
We obtain similar behavior and significant gradients near the walls. Additionally, our DDFV scheme takes approximately 30 minutes to solve on a laptop (CPU: Intel(R) Xeon(R) Gold 6254 CPU @ 3.10GHz; MEMORY: DIMM DDR4 Synchronous Registered (Buffered) 2933 MHz (0.3 ns)), compared to 1 hour in [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF]. (red). Black curves represent linear trends of the potential V in the bulk.

We observe that outside of the boundary layer, the ionic concentration c P and c N are very close, which corresponds to an electroneutral bulk. We compute the absolute value of this difference at time t = 0.5: q K * ∈M * := (|c K * ,50 P c K * ,50 N |) K * ∈M * . In [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF], the authors reported a maximum value for this difference of 3.6 × 10 -6 in the interval [0.25, 0.75] at t = 0.5. We plot the values for q K * ∈M * cut at 3.6 × 10 -6 , and color the regions where q > 3.6 × 10 -6 in black (Fig. 6A). We observe that in our simulations, the x-boundary of the region such that q K * ∈M * ≤ 3.6 × 10 -6 is [0.075, 0.925], which is a larger region than the segment [0.25, 0.75] in [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF]. Note that the maximum for q K * ∈Z * [0.25,0.75]×{0.5} on [0.25, 0.75] × {0.5} is 1.5 × 10 -6 , which is half the value 3.6 × 10 -6 observed in [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF]. In Fig. 6B, we plot the values of q K * ∈Z * [0.075,0.925]×{0.5} (see ( 10)). The error is decreasing from the boundaries (x = 0 and x = 1) to the center x = 0.5. At x = 0.5, the error is reaching the machine error 10 -14 . We show with this test case that our DDFV scheme can capture well the boundary layer dynamics. Compared to [START_REF] Song | Electroneutral models for dynamic poisson-nernst-planck systems[END_REF], our electroneutral zone is larger, suggesting that the perturbation of the numerical results due to the high gradients inside the boundary layer, is less pronounced with our scheme.

Applications to neuroscience

In the last section, we apply our DDFV scheme to investigate the dynamics of voltage and ionic concentration in different neuronal geometries. We consider two specific geometries in 2D: a bifurcation in the dendritic tree, and a dendritic spine. The dendritic tree has a treelike geometry, and a bifurcation corresponds to the location where one branch divides into two branches (Fig. 7). Dendritic spines are mushroom-like protrusions of a few micrometers, found on the dendritic tree. In section 5.1, we consider a dendrite bifurcation with two thin branches connecting to a larger one and simulate the propagation of an influx of ions arriving at the edge of one of the thin branches. We investigate three scenarios: the influence of an ionic reservoir on voltage and ionic concentration dynamics, the summation of two simultaneous influx of ions and the propagation of a single influx in the rest of the domain. Section 5.2 focuses on dendritic spines. We first investigate the influence of the spine head geometry on voltage and ionic concentration dynamics. We then consider two neighboring spines and investigate the invasion of an influx of ions arriving in one spine on the other.

We recall here the definitions of the initial and boundary conditions given in section 2 (values are given in Table 3):

c P (t = 0, x) = c 0 P , c N (t = 0, x) = c 0 N , V (t = 0, x) = V 0 ∀x ∈ Ω. ( 30
)
We impose constant Dirichlet boundary condition on Γ Dir :

c P = c 0 P , c N = c 0 N , V = V 0 , on Γ Dir × (0, T f ).
We consider ∂Ω r for homogenous and ∂Ω i for non homogeneous Neumann boundary conditions, such that Γ N eu = ∂Ω r ∪ ∂Ω i . The homogeneous Neumann boundary condition on ∂Ω r models the impermeability of the neuronal membrane to ions. The non-homogeneous boundary condition on ∂Ω i models the influx of ions I received by the dendrite, at specific locations called synapses:

∇V • n = 0, on ∂Ω i ∪ ∂Ω r ×(0, T f ), ∇c N • n = 0, on ∂Ω i ∪ ∂Ω r ×(0, T f ), ∇c P • n = 0, on ∂Ω r ×(0, T f ), ∇c P • n = I, on ∂Ω i ×(0, T f ).
(31)

I(t) = I stim (t) πr 2 i F D P
, where

I stim (t) = I max t τ exp - t τ + 1 (32)
represents the injected synaptic current, inspired by [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF].

Note that in this section, the source terms f cp , f cn and f V are set to zero. We define a time t 0 for each simulation, corresponding to the time at which the concentration c P is maximal on Ω. Note that it corresponds to a maximum in both space and time. We precise the value of t 0 in each case. We use electrodiffusion parameters as defined in Table 3. Finally, we present in this section the numerical results for c P , as the results for c N are qualitatively similar. [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF] Table 3: Electrodiffusion parameters.

F

Propagation of an influx of ions at a dendritic tree bifurcation

In this part, we model signal propagation at a bifurcation in the dendritic tree. We consider a domain Ω B representing the bifurcation, where a large branch divides into two thinner ones (Fig. 7). We model the influx of ions as a current I(t) (eq. ( 32)) injected at the end of the two thin branches (∂Ω up i for the upper branch, and ∂Ω down i for the lower branch, Fig. 7,

∂Ω i = ∂Ω up i ∪∂Ω down i
). We impose Dirichlet boundary conditions at the end of the large branch (Γ Dir on Fig. 7), to represent the connection with a larger dendrite. This large dendrite is considered an ionic reservoir due to its large size, i.e. it has fixed ionic concentrations. The different lengths and nodes defining the domain Ω B are described in Fig. 7, with the length values and node coordinates given in Table 4.

We define two positions near the end of each branch, namely (I) and (J), to compare the voltage and concentration dynamics in the branches. We use a triangular mesh with 6655 cells and a mesh size h = 0.36. The simulations are performed with a final time of T f = 1.5 sec, and a time step of dt = 5 × 10 -3 sec. 4: Geometric parameters for domain Ω B , representing a dendritic bifurcation.

Effect of the distance to an ionic reservoir on voltage and concentration dynamics

Our model considers that the end of the large branch is connected with a wider one, called an ionic reservoir, such that the ionic concentrations at Γ Dir are constant. The distance between this ionic reservoir and the bifurcation influences the dynamics of voltage and ionic concentrations everywhere on the domain Ω B . We expect that the longer the large branch (i.e. the longer L 5 , Fig. 7), the smaller the influence on the dynamics close to the bifurcation point. To measure this influence, we apply a current at the end of the two thin branches (∂Ω i on Fig. 7).

We then realize simulations for different values of L 5 ranging from 11 µm to 33 µm (see Fig. 8 for the case L 5 = 33 µm). We finally compute the difference between each numerical solution, and consider this difference relative to the peak amplitude. We consider that the influence of the Dirichlet boundary condition on ionic and voltage dynamics is small when this relative difference is below 3.6%. We realize numerical simulations in domain Ω B for three different values of L 5 : 11 µm, 22 µm and 33 µm. Note that for each value of L 5 , the domain is modified, and so is the mesh. We indicate the different coordinates and parameter values that are modified in the three configurations of domain Ω B in Table 5. To account for the mesh modifications in the discrete solution for the different configurations, we note c T ,n P,q the discrete value of the concentration c P in the case L 5 = q µm. We compare the numerical solutions obtained on different domains, by evaluating their difference on two straight lines y 1 : [0, 11] × {1} (see Fig. 7), and y 2 : [-11, 11] × {1}. Note that the line y 2 extends the line y 1 for the configurations L 5 = 22 µm and L 5 = 33 µm. We note e K * ,n p,q K * ∈Z * y i the absolute value of the difference of the dual values for c P , on line y i (see ( 10)), for the domains with L 5 = p and L 5 = q, at time iteration n, for i = 1, 2:

e K * ,n p,q K * ∈Z * y i = c K * ,n P,p K * ∈Z * y i -c K * ,n P,q K * ∈Z * y i . ( 33
)
Finally, we compute the percentage difference relative to A, between the simulation results on line y i for t = 0 to t = 0.5 sec (time steps n = 0, .., 100), and plot the result in We observe that the difference between the solutions is always maximal at the Dirichlet boundary condition, at time t = 0, (Fig. 9). The maximal percentage difference between the concentration values c P for L 5 = 11 µm and L 5 = 22 µm is 9.1 % (A). For L 5 = 11 µm and L 5 = 33 µm, the maximum percentage difference is 10.4 % (B), and finally, for L 5 = 22 µm and L 5 = 33 µm, to 3.6 % (C). Finally, we conclude that as expected, the influence of the Dirichlet boundary condition Γ Dir (i.e. the ionic reservoir) on concentration dynamics c P , decreases when the length L 5 of the large branch increases. Specifically, for our geometry, the influence is small if the reservoir is at a distance larger than 22 microns. -11, 11], so as the x-axis in panel C.

Signal summation

We now study signal summation at a bifurcation in the dendritic tree, when the ionic reservoir is far from the bifurcation. We compare ionic concentration and voltage dynamics in two different scenarios: when only the upper branch receives the influx (∂Ω up i in Fig. 7, scenario 1, ∂Ω down i is turned into a homogeneous Neumann boundary condition) and when there is an influx of ions at both ends of the two thin branches (∂Ω up i and ∂Ω down i in Fig. 7, scenario 2). We will call a branch 'active' if it receives an influx of ions, and 'inactive' if it does not. We run simulations in the domain with L 5 = 33 µm (Fig. 7 andTable 5). We consider the line y 3 = [-22, 11] × {1}, and the set of dual cells Z * y 3 (see ( 10)). In Fig. 10, we compare the concentration (c K * ,21

P -c 0 P ) i K * ∈Z * y 3 and voltage (V K * ,21 ) i K * ∈Z * y 3
at time t 0 = 0.105 sec, for the two scenarios 1 and 2. In both graphs, the brown curve represents scenario 1 with one active branch, while the blue curve corresponds to scenario 2 where the two branches are active. We observe that, on line y 3 , the concentration and potential values in scenario 2 (blue) are twice the values in scenario 1 (brown). , at peak time t 0 = 0.105 sec, for i = 1 in brown and i = 2 in blue. In dashed red, we plot 2(c K * ,21

P -c 0 P ) 1 K * ∈Z * y 3 . B: Dual values (V K * ,21 ) i K * ∈Z * y 3
, at peak time t 0 = 0.105 s, for i = 1 in brown and i = 2 in blue. In dashed red, we plot 2(V K * ,21 ) 1

K * ∈Z * y 3
.

In Fig. 11, we plot the time evolution of the concentration c Pc 0 P and the potential V at position (H), i.e. (c

K * H ,n P -c 0 P ) i n=0,••• ,N and (V K * H ,n ) i n=0,••• ,N
for the two scenarios. We recall that the definition of K * H is given by (11). As in Fig. 10, we observe that the potential and concentration dynamics in scenario 2 are twice the dynamics in scenario 1.

K * H ,n P -c 0 P ) i n=0,••• ,N for i = 1 in brown and i = 2 in blue. In dashed red, we plot 2(c K * H ,n P -c 0 P ) 1 n=0,••• ,N max n=0,••• ,N (c K * H ,n P ) 2 -2(c K * H ,n P ) 1 = 5 × 10 -3 . B: Dual values (V K * H ,n ) i n=0,••• ,N for i = 1 in brown and i = 2 in blue.
In dashed red, we plot

2(V K * H ,n ) 1 n=0,••• ,N . We have max n=0,••• ,N (V K * H ,n ) 2 -2(V K * H ,n ) 1 = 0.3 .
We finally compute the absolute difference between the concentration (c

K * ,n P -c 0 P ) 2 K * ∈Z * y 3 n=0,••• ,N
and twice the concentration

(c K * ,n P -c 0 P ) 1 K * ∈Z * y 3 n=0,••• ,N
on line y 3 and for t ∈ [0, T f ]. This difference is zero at the Dirichlet boundary condition (x = -22 µm, for all t), and stays below 10 -2 mM up to x = 8 µm (resp. 10 -1 mV up to x = 8 µm) for the concentration (resp. the potential). It then steeply increases to 0.33 mM (resp. 0.35 mV) for the concentration (resp. the potential) close to x = 11 µm (bifurcation point (H)). In summary, the signal in the large branch is doubled when two branches are active compared to only one active branch. This result is coherent with a linear summation of the signal in passive dendrites observed experimentally [START_REF] Cash | Linear summation of excitatory inputs by ca1 pyramidal neurons[END_REF].

Branch invasion

We now investigate the invasion of a signal in an inactive branch, by plotting the evolution of the concentration and potential dynamics in the small branches in scenario 1 (only the upper branch is active). We consider the nodes (I) and (J) toward the end of the upper and lower branches respectively, far from the bifurcation point (See Fig. 7 and Table 5). In Fig. 12, we plot the time evolution of the concentrations (c

K * I ,n P -c 0 P ) 1 n=0,••• ,N and (c K * J ,n P -c 0 P ) 1 n=0,••• ,N and the potential (V K * I ,n ) 1 n=0,••• ,N and (V K * J ,n ) 1 n=0,••• ,N
at the points (I) and (J) (see eq. ( 11)). We observe that the influx of ions arriving at the upper branch creates a transient rise in the concentration and in the potential that invades the lower branch. We compare the time and amplitude of c P and the potential at the two positions (I) and (J) and observe that the maximum of the signal arrives in (J) with a delay of 0.21 sec, and an amplitude reduced by 82 %, decreasing from 81.43 mM in node (I) to 14.57 mM in node (J) (A). Interestingly, the concentration curves overlap at (I) and (J), for t > 0.5 sec. We hypothesize that this is due to the vanishing of the potential at t > 0.5 sec, highlighting the role of the convective term in ionic dynamics (ion transport resulting from the electric field). Concerning the potential dynamics, the maximum of the signal arrives in (J) with no delay and a decrease in amplitude of 31 %, decreasing from 44.81 mV to 30.91 mV. In summary, we observe a discrepancy between a rapid and strong invasion of the electrical signal in the inactive branch, and a low invasion of the ionic concentration signal, with a delay of a few hundred ms. A: (c

K * J ,n P -c 0 P ) 1 n=0,.
.,N , at position (J) in red and (c

K * I ,n P -c 0 P ) 1 n=0,..,N , at position (I) in blue. B: (V K * J ,n ) 1 n=0,.
.,N , at position (J) in red and (V K * I ,n ) 1 n=0,..,N , at position (I) in blue.

Modeling and simulation of ionic and voltage dynamics in dendritic spines

In this section, we investigate the effect of the specific geometry of dendritic on voltage and ionic dynamics. Dendritic spines are located on the dendritic tree and serve as sites for receiving synaptic input in the form of an influx of ions. Their variations in size and shape are usually associated with the neuronal coding of learning and memory. Extensive efforts have been made to develop experimental techniques for visualizing and analyzing ionic and electric field dynamics in dendritic spines. However, due to their microscopic scale, such experimental investigations remain challenging.

Influence of the spine head geometry on ionic and voltage dynamics

Dendritic spines have a bulbous head connected to a thin neck, that we model in 2D by a domain Ω S composed of a circle with radius r (the head) connected to a rectangle of length L 1 and width D 1 (the neck) (Fig. 13). All parameters related to domain Ω S are defined in Table 6. The influx of ions is modeled, as in the previous section, with a non-homogeneous Neumann boundary condition on ∂Ω i . We set, as previously described, homogeneous boundary conditions on ∂Ω r , as well as Dirichlet boundary conditions on Γ Dir (see Fig. 13).

The mesh used in our simulations is gradually refined along the boundary ∂Ω S of domain Ω S and consists of 6589 triangular cells with a mesh size of h = 7.33 × 10 -2 . As an example, Figure 13 shows a zoomed-in view of the mesh at the junction between the head and the neck. 0.2 µm Diameter of the neck r i 0.04 µm Length of ∂Ω i (Fig. 13) Table 6: Geometric parameters for domain Ω S , representing a dendritic spine

In the following, we present the numerical simulations of the scheme ( 15)-( 25) on domain Ω S , with boundary conditions (31). The simulation is performed with a final time of T f = 0.5 sec and a time step of dt = 5 × 10 -3 sec. In Fig. 14, we plot the ionic concentration and voltage dynamics at time t 0 = 0.075 sec in domain Ω S . Panel A (resp. panel B) shows the dual concentration values (c K * ,15

P -c 0 P ) K * ∈M * (resp. (V K * ,15 ) K * ∈M * of the potential V ), at time t 0 = 0.075 sec.
Figure 15: Schematic representation of a dendritic spine, illustrating the lines {x = 0} (dashed black) and {y = k} (dashed red) for k ∈ {0.5, 1, 1.1, 1.3, 1.5, 1.7, 1.9, 2}. The points P k are the intersections between the lines {x = 0} and {y = k} (red star).

We plot in Fig. 16 the time evolution of the dual concentration value (c

K * P k
,n P c 0 P ) n=0,..,N (A) and the dual potential value (V

K * P k
,n ) n=0,..,N (B) at points P k (eq. ( 11)). We observe that the peak amplitude for c Pc 0 P (resp. V ) is equal to 38.59 mM (resp. 6.37 mV) at (P 2 ) and decreases to 26.16 mM (resp. 4.36 mV) at (P 1 ), i.e. a decrease of more than 30 % for both the concentration and voltage. We also observe that the decrease is faster in the neck, with the peak amplitude reaching 13.35 mM (resp. 2.31 mV) at P 0.5 for the concentration (resp. voltage). In [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF] the peak amplitude of c P -c 0 P (resp. V ) at the point representing the head, reaches around 33 mM (resp. around 5 mV) (see Fig. 4C, resp. Fig. 2D-E in [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF]), which is very close to our 2D results at point P 2 . Nevertheless, our results indicate that almost a third of the signal is lost within the head, highlighting the necessity to simulate the dynamics inside the head. ,n ) n=0,..,N .

We finally investigate the spatial variation of c Pc 0 P within the spine, at peak time t 0 =0.075 sec. In Fig. 17A, we plot the dual concentration values (c K * ,15 P c 0 P ) K * ∈Z * x=0 (see eq. ( 10)). We observe that the concentration decreases linearly within the neck, reaching the Dirichlet boundary condition at (x, y) = (0, 0), which suggests that a 1D approximation is relevant in the neck. Nevertheless, in the head, the variation of the concentration profile is more complex, with a higher gradient at the top of the head, near the influx of ions (∂Ω i ), followed by a lower gradient down to the neck. In Fig. 17B, we plot (c K * ,15 P c 0 P ) K * ∈Z * y=k for k ∈ {0.5, 1, 1.3, 1.5, 1.7, 1.9, 2.0} (eq. ( 10)). The variations along the lines Y 0.5 and Y 1 , located inside the neck, are very small. In contrast, the geometry of the head creates concentration gradients, that are more pronounced close to the boundary receiving the influx of ions, and close to the head-neck junction. On line Y 1.9 , the concentration decreases from 35.46 mM at position x = 0, to 33.33 mM at the boundary (black curve on Fig. 17D), which corresponds to a decrease of 6% relative to the reference value c In summary, we observe that the ionic and voltage dynamics in the neck of dendritic spines are well captured by a 1D model. In contrast, to capture the specific dynamics of voltage and concentration inside the head, one needs to take into account its bulbous geometry. In particular, in our 2D simulations, we observe that the peak amplitude is decreased by almost 30 % between the top of the head and the head-neck junction, and the decrease is non-linear. This justifies the need for at least a 2D approach to obtain accurate simulations of ionic concentration dynamics in a spine.

5.2.2

Influence of the distance between an ionic reservoir and a dendritic spine on the ionic and voltage dynamics.

In this part, we consider two neighboring dendritic spines and investigate the influence of a signal arriving in one spine, on the ionic and voltage dynamics of the other (Fig. 18). We test the hypothesis that dendritic spines act as autonomous compartments, isolating the material located at their head from the rest of the dendritic tree [START_REF] Hering | Dentritic spines : structure, dynamics and regulation[END_REF], versus the hypothesis of signal invasion, when the voltage and concentration in a spine is substantially increased by an ionic influx arriving in a neighboring one. We consider the domain Ω 2S (Fig. 18), with parameters given in Table 7. Within one of the spines, which we call the active spine, we model an influx of ions at boundary ∂Ω i , using a non-homogeneous Neumann boundary condition, similar to the previous subsections, eq. ( 32). The spine that does not receive an influx of ions is called the inactive spine. We address the questions of the distance at which the influx of ions propagates from the active spine, and to what extent the inactive spine perceives this influx. We use a triangular mesh with 2385 cells and a mesh size h = 0.16. The simulations are performed with a final time of T f = 5 sec, and a time step of dt = 5 × 10 -3 sec.

(A) (0,1.5) Coordinates of point (A) (B) (0,0) Coordinates of point (B) (C) (1,-0.25) Coordinates of point (C) (D) (2,0) Coordinates of point (D) (E) (2,1.5) Coordinates of point (E)
Table 7: Parameters for the domain Ω 2S , representing two neighboring dendritic spines (Fig. 18).

We investigate the effect of the distance L 3 , between a spine and the Dirichlet boundary condition Γ Dir , on voltage and ionic dynamics. We consider several domains with L 3 ranging from 0.9 µm to 23.9 µm. In each configuration, the length L 4 is modified such that L 4 = 2L 3 + 2D 1 + L 2 , and L 2 is set to 1.8 µm (Fig. 18). The different mesh information are given in Table 8.

Value for L 3 0.9 µm 2.9 µm 5.9 µm 11.9 µm 23.9 µm Fig. 19 illustrates the influence of the length L 3 on the dynamics of c Pc 0 . We plot the dual values (c K * ,20 P c 0 ) K * ∈M * , (c K * ,20 P c 0 ) K * ∈M * and (c K * ,20 P c 0 ) K * ∈M * at times t 0 = 0.1 sec, for three domains Ω 2S with L 3 = 0.9 µm, 2.9 µm and 5.9 µm. Note that the minimum value for c Pc 0 is imposed by the Dirichlet boundary condition at Γ Dir (c Pc 0 = 0 mM, with c 0 = 163 mM, Table 3). We observe that the size of the domain influences the values of the concentration within the two heads, with a larger L 3 leading to a higher concentration in the active spine, and a lower concentration value in the inactive one. We also observe a slight increase in c Pc 0 peak time when L 3 increases. Hence, the results suggest that the farther the spine is from a large ionic reservoir, the farther the signal propagates. We plot the dual values (c K * ,20 P c 0 ) K * ∈M * in domain Ω 2S , with L 3 = 0.9 µm (A), L 3 = 2.9 µm (B) and L 3 = 5.9 µm (C) at t 0 = 0.1 sec. The range of the colorbar is fixed between plots, with a maximum value c Pc 0 = 180.2 mM, corresponding to the maximum of c Pc 0 in the configuration with L 3 = 5.9 µm. We then investigate the propagation of the concentration within the dendrite and inactive spine at several specific points (A) -(E) (Fig. 18, Table 7) in the five different configurations (L 3 ranging from 0.9 µm to 23.9 µm). We plot in Fig. 20 the time evolution of the concentrations (c K * X ,n P c 0 P ) n=0,..,N (see ( 11)). Note that the curves plotted in D and E almost superimpose, with a difference between them lower than 1 mM. We observe the effect of the Dirichlet boundary condition on the time evolution of the concentration. As expected, the shorter L 3 , the bigger the impact on signal propagation. Indeed, for L 3 = 0.9 µm (red), we observe that the max of c Pc 0 goes from 109.62 mM in (A) to 6.60 mM in (B), and is below 10 -5 mM at (C), (D) and (E), i.e. the signal does not propagate inside the dendrite and the inactive spine. This means that for short L 3 the active spine behaves as an autonomous compartment, isolating its material from the rest of the dendritic tree. For L 3 = 2.9 µm, the max of c Pc 0 goes from 145.95 mM in (A) to 55.92 mM in (B), and 24.90 mM in (E). Hence, the ratio of signal reaching (E), defined as r c

A

-E = max n=0,..,N c K * E ,n P -c 0 P max n=0,..,N c K * A ,n P -c 0 P , is 17
%. For L 3 = 11.9 µm and 23.9 µm, a threshold seems to emerge, with 32 % of the peak concentration in (A) transmitted to (E). We also observe a shift in the peak time between position (A) and (E), denoted by ∆ c

A-E , going from 0.03 sec for L 3 = 0.9 µm, to 0.09 sec for L 3 = 2.9 µm and 0.16 sec for L 3 = 23.9 µm. We observe similar behavior for the voltage, where the ratio of the signal reaching (E), r

V A-E = max n=0,..,N V K * E ,n max n=0,..,N V K * A ,n , is 8% (resp. 63%, resp .
81%) for L 3 = 0.9 µm (resp. L 3 =11.9 µm, resp. L 3 =23.9 µm). The shift in the peak time between the electrical signal in (A) and the electrical signal in (E), denoted by ∆ V A-E , is below 0.02 sec in all curves. Table 9 gives, for each configuration, the time shift ∆ c

A-E , the peak amplitude max A-E of signal reaching (E), for the concentration dynamics, as well as the time shift ∆ V A-E , the peak amplitude max n=0,..,N V K * E ,n at position (E) and the ratio r V A-E of signal reaching (E) for the potential dynamics. This indicates that in this condition, a signal arriving in a spine is influencing the ionic concentration and voltage in its inactive neighbors, which we call a signal invasion. The voltage invasion is more important than the ionic concentration invasion. To summarize, our simulations suggest that depending on the distance to the closest ionic reservoir, a spine can either act as an autonomous compartment isolated from its neighbors or be subject to signal invasion. V K * E ,n at position (E) and ratio r V A-E of the electrical signal from (A) reaching (E), for the five different configurations with L 3 going from 0.9 µm to 23.9 µm.

L 3 ∆ c A-E max n=0,..,N c K * E ,n P -c 0 P r c A-E ∆ V A-E max n=0,..,N V K * E ,n r V A-E 0.
In Table 11 we give, for each configuration, the time shift ∆ c

A-E , the peak amplitude max A-E of signal reaching (E) for the concentration dynamics, as well as the time shift ∆ V A-E , the peak amplitude max n=0,..,N

V K * E ,n at position (E) and ratio r V A-E of signal reaching (E) for the potential dynamics. These results suggest that spines receive electrical signals from other spines far away, as a spine at a 10.8 µm distance still receives 63 % of the signal, with tens of ms of delay. The variation in concentration is sensed to a lesser extent, as a spine at a 10.8 µm distance only perceives 14% of the variation of concentration, with a delay superior to 0.5 sec. V K * E ,n at position (E) and ratio r V A-E of the electrical signal from (A) reaching (E), for the five different configurations with L 2 going from 0.9 to 10.8 µm.

Discussion and conclusion

In this paper, we present an algorithm for simulating the Poisson-Nernst Planck system of equations in two dimensions using the DDFV method, while preserving the positivity of ionic concentration. We challenge our algorithm using different test cases and achieve a secondorder accuracy in space, which is consistent with existing results coming from literature. We also show that our system can capture a boundary layer when present. The proof of the existence of the discrete solution of ( 15)-( 25) is both too long and not the goal of the present paper. Hence, a following paper will detail the existence of a solution, as well as positivity and convergence, most probably using similar techniques as in [START_REF] Moatti | A structure preserving hybrid finite volume scheme for semiconductor models with magnetic field on general meshes[END_REF]. Using our DDFV framework, we then investigate the ionic and voltage dynamics in twodimensional specific geometries of the neuronal dendritic tree: a branch bifurcation and a dendritic spine. First, our two-dimensional results show that the voltage and ionic dynamics in a dendritic spine are not well approximated by one-dimensional models. Indeed, if a one-dimensional approximation for both voltage and concentration is relevant in the neck, the bulbous geometry of the head requires at least a two-dimensional approach, as the space dynamics is highly non-linear, with steep gradients at the top of the head and around the head-neck junction. Our numerical results also highlight that dendritic spines can sense electrical signals far away on a thin branch, whereas it is not the case for ionic concentration transients, that reach only close spines. In our configuration, a spine at a distance 10.8 µm of a spine receiving a signal will get 62% of the electrical signal and only 14% of the concentration signal. One direct consequence of this, is that calcium imaging cannot be used as a tool to observe sub-threshold electrical signal propagation in neuronal compartments, especially when voltage-gated calcium channels are not opening. The actual development of voltage sensors could resolve this diffi-culty. Our final result concerns the influence of an ionic reservoir on voltage and calcium dynamics. This reservoir can represent any large branch that the neuronal compartment is connected to. The rationale behind this is that the large branch being so large compared to the thin branch, we can consider the ionic concentrations and voltage to be constant at the connection. In the equations, the reservoir is modeled using a Dirichlet boundary condition and behaves as a sink. Our simulations show that the close proximity of an ionic reservoir such as the dendritic shaft or any large compartment, is killing the signal, preventing it from propagating and invading close neuronal structures such as thin branches at a bifurcation or dendritic spines. On the opposite, a signal arriving at the leading edge of the dendritic tree, far from an ionic reservoir and where only small branches are present, will propagate at a larger distance and invade neighboring dendritic spines. From the spine point of view, these results show that the same spine at different positions in the tree would behave differently: a spine close to the dendritic shaft would more likely act as an autonomous compartment, compared to a collection of spines located in small dendritic protrusions, that would be keener to signal invasion and to influence each other. Hence, in addition to the geometry of the spine, the local geometry of the dendritic tree is shaping spine behavior, in that the same spine at different positions in the tree would behave differently. Hence, the position of a dendritic spine relative to the entire tree (i.e. close to the soma versus at the distal edge of the tree) shapes its function, making plasticity not at the level of the spine, but at the level of the full dendritic geometry. These results are in line with several experimental observations showing that synaptic development is spatially regulated inside the dendritic tree [Druckmann et al., 2014, London and[START_REF] London | Dendritic computation[END_REF]. In particular, the electrical compartmentalization of dendritic spines is nowadays a leading question in developmental neuroscience, with various experimental studies showing a variety of results depending on the types of cells or on their developmental stages, and going from signals spreading locally and invading neighboring spines to compartmentalization [START_REF] Vh Cornejo | Voltage compartmentalization in dendritic spines in vivo[END_REF][START_REF] Yuste | Dendritic spines[END_REF][START_REF] Kfh Lee | Correlated synaptic inputs drive dendritic calcium amplification and cooperative plasticity during clustered synapse development[END_REF][START_REF] Yuste | Dendritic spines and distributed circuits[END_REF]. Our study would suggest considering these differences in relation to the actual size of the dendritic tree, and especially the presence of large dendrites. In the same way, the distance to the nearest large dendritic compartment should be considered while investigating cooperative, homosynaptic and heterosynaptic plasticity [START_REF] Te Chater | My neighbour hetero-deconstructing the mechanisms underlying heterosynaptic plasticity[END_REF]Goda, 2021, Sjöström and[START_REF] Sjöström | A cooperative switch determines the sign of synaptic plasticity in distal dendrites of neocortical pyramidal neurons[END_REF].

Figure 2 :

 2 Figure 2: Left: Diamond mesh D (green). One cell is colored in light green. Right: A diamond D (green) with direct basis (τ K * ,L * , n σK ) and (n σ * K * , τ K,L ). σ (resp. σ * ) is the primal (resp. dual) edge.

Finally

Figure 3 :

 3 Figure 3: Example of a quadrangle mesh with h = 2.341 × 10 -1 .

Figure 4 :

 4 Figure 4: Simulation results of the third test case adapted from [Song et al., 2018] at T f = 1. A: (c K * ,100 P

Figure 5 :

 5 Figure 5: A: Dual discrete values of c P and c N on line y = 0.5 at time t = 0.5, (c K * ,50 P

Figure 6 :

 6 Figure 6: A: Absolute difference between c P and c N at time t = 0.5, q K * ∈M * = (|c K * ,50 P c K * ,50 N |) K * ∈M * . q K * ∈M * goes from 10 -14 to 3.6 × 10 -6 . Regions where q K * ∈M * > 3.6 × 10 -6 are colored in black. B: Plot of q K * ∈Z * [0.25,0.75]×{0.5} on the segment [0.075, 0.925] × {0.5}, in logarithmic scale.

Figure 7 :

 7 Figure 7: Geometry of the bifurcation domain Ω B , in the configuration L 5 = 11 µm. The specific line y 1 = [0, 11] × {1}, where we monitor concentration and voltage dynamics is plotted in dashed red. The coordinates of each node are given in Table4.

  × 10 -1 3.9 × 10 -1 4 × 10 -1 Position (x,y) of node (DDifferent coordinates and parameter values that are modified in the three configurations of domain Ω B : L 5 = 11 µm, L 5 = 22 µm and L 5 = 33 µm.In Figure8, we plot the dual values c K * ,21 P,33c 0 P K * ∈M * , corresponding to the concentration c Pc 0 P at time t 0 = 0.105 sec, in the case L 5 = 33 µm. We observe that the maximum of the solution, A = maxK * ∈M * c K * ,21P,33c 0 P = 93.8 mM, is reached at the injection boundary ∂Ω i . Along x, the solution decreases to reach the value c 0 P at the Dirichlet boundary Γ Dir .

Figure 8 :

 8 Figure 8: Dual values c K * ,21 P,33c 0 P K * ∈M * of the concentration c Pc 0 P at time t 0 = 0.105 sec, in the case L 5 = 33 µm.

Fig

  Fig. 9: e 1 :=

Figure 9 :

 9 Figure 9: Comparison of the numerical solution between the three domains with L 5 = 11 µm, L 5 = 22 µm and L 5 = 33 µm. A (resp. B, resp. C): Plots of e 1 (resp. e 2 , resp. e 3 ) the absolute values of the differences of the dual values for c P , on line y 1 (resp. y 1 , resp. y 2 ), for the domains with L 5 = 11 µm and L 5 = 22 µm, (resp. L 5 = 11 µm and L 5 = 33 µm, resp. L 5 = 22 µm and L 5 = 33 µm). Note that line y 2 ranges in [-11, 11], so as the x-axis in panel C.

Figure 10 :

 10 Figure 10: Evolution of the concentration and voltage dynamics on line y 3 for the two scenarios 1 and (resp. 2), with one (resp. two) active branch(es). A: Dual values (c K * ,21 P

Figure 11 :

 11 Figure 11: Evolution of the concentration and voltage dynamics at point (H) for the two scenarios 1 and 2. A: Dual values (c

Figure 12 :

 12 Figure 12: Time evolution of the dynamics of c P and V in the small branches for scenario 1. A: (c

Figure 13 :

 13 Figure 13: A: Domain Ω S representing a dendritic spine. B: Mesh for domain Ω S zoomedin at the head-neck junction where it is refined close to the boundaries. The mesh size is h = 7.33 × 10 -2 .

Figure 16 :

 16 Figure16: Time evolution of the ionic concentration c Pc 0 P and voltage V in the dendritic spine, at position P k , for k ∈ {0.5, 1, 1.3, 1.5, 1.7, 1.9, 2}. A: (cK * P k ,n

Figure 17 :

 17 Figure17: A: Spatial evolution of (c K * ,15 P c 0 P ) K * ∈Z * x=0 on line x = 0 at time t 0 = 0.075 sec. B: Spatial evolution of (c K * ,15 P c 0 P ) K * ∈Z * y=k for k ∈ {0.5, 1, 1.3, 1.5, 1.7, 1.9, 2.0} at time t 0 = 0.075 sec.

Figure 18 :

 18 Figure 18: Domain Ω 2S , representing two neighboring dendritic spines. Each spine has a head (circle with centers (A) and (E)) and a neck (thin vertical rectangle). They are both connected to a dendrite (large horizontal rectangle, with points (B), (C) and (D)).

Figure 19 :

 19 Figure 19: Spatial variations of the concentration c P in three domains with increasing L 3 . We plot the dual values (c K * ,20

  Figure 20: Time evolution of (c K * X ,n P c 0 P ) n=0,..,N , for different values of L 3 ranging from 0.9 to 23.9 µm, at positions (A)-(B)-(C)-(D)-(E). A: (c K * A ,n P c 0 P ) n=0,..,N . B: (c K * B ,n P c 0 P ) n=0,..,N . C: (c K * C ,n P c 0 P ) n=0,..,N . D: (c K * D ,n P c 0 P ) n=0,..,N . E: (c K * E ,n P c 0 P ) n=0,..,N .

P

  at position (E) and the ratio r c

  position (E), ratio r c A-E of the concentration signal in (A) reaching the inactive spine (E), time shift ∆ VA-E of the electrical signal, peak amplitude max n=0,..,N

P

  at position (E), and the ratio r c

P

  of the concentration at position (E), ratio r cA-E of the concentration signal in (A) reaching the inactive spine (E), time shift ∆ VA-E of the electrical signal, peak amplitude max n=0,..,N

Table 1 :

 1

	4 1.67	4.713E-4 1.61	5.137E-3 2.00
	1.212E-1 6.25E -4 2.066E-4 2.08 6.253E-2 1.56E -4 4.897E-5 2.18	1.220E-4 2.05 2.843E-5 2.20	1.312E-3 2.07 3.182E-4 2.14

Table 8 :

 8 Parameter values that are modified in the five configurations of domain Ω 2S .

Table 9 :

 9 Time shift ∆ c A-E in the concentration signal, peak amplitude max

	9 µm	0.03 sec	10 -5 mM	0 % 0.01 sec	1.65 mV	8 %
	2.9 µm	0.09 sec	24.9 mM	17 % 0.01 sec	7.62 mV	26 %
	5.9 µm	0.09 sec	40.94 mM	27 % 0.02 sec	17.58 mV	44 %
	11.9 µm 0.16 sec	52.44 mM	32 % 0.02 sec	38.41 mV	63 %
	23.9 µm 0.16 sec	53.45 mM	32 % 0.02 sec	103.46 mV	81 %

Table 11 :

 11 Time shift ∆ c A-E in the concentration signal, peak amplitude max

	r V A-E

We now investigate the effects of the specific geometry of the spine on ionic concentration and voltage dynamics, and in particular the spatial variations of voltage and ionic concentration in the spine head. We then compare our 2D-simulation results with the 1D results presented in [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF]. In [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF], the authors introduced a temporal deconvolution procedure (STAR method) to recover voltage dynamics in dendritic spines, from fluorescence images of a genetically-encoded voltage sensor. This procedure is based on the PNP system of equations. From a domain similar to Ω S , the authors derived a 1D equivalent model reducing the spine head as a point and simulated the PNP system of equations using COMSOL. Note that we use the same input current at ∂Ω i as in [START_REF] Cartailler | Deconvolution of voltage sensor time series and electro-diffusion modeling reveal the role of spine geometry in controlling synaptic strength[END_REF].

We consider several lines on the domain Ω S : {x = 0} and the lines {y = k} for k ∈ {0.5, 1, 1.3, 1.5, 1.7, 1.9, 2} (Fig. 15). We finally define the points P k , at the intersection between {x = 0} and {y = k}, i.e. the points (0, k).

Influence of the distance between two spines on signal invasion

This part focuses on varying the distance L 2 between two spines (Fig. 18). We consider five different values of L 2 , ranging from 0.9 to 10.8 µm, and set L 3 to 23.9 µm, to ensure that we reach the plateau in signal invasion observed in subsection 5.2.2, for all the configurations. The mesh information are given in Table 10 We observe, as expected, that concerning the variation of concentration, the closer the two spines, the larger the signal received by the inactive one. The shift in the peak time between position (A) and (E), ∆ c

A-E , increases when the distance between the spines L 2 increases. We observe a significant shift in the peak concentration time in the inactive spine compared to the active one, going from ∆ c

A-E = 0.145 sec in the case of close spines (L 2 = 0.9 µm) to ∆ c

A-E = 0.535 sec for spines that are farther away (L 2 = 10.8 µm). The ratio r c A-E of signal reaching (E) decreases while increasing L 2 , going from 37 % when L 2 = 0.9 µm, to 14 % for L 2 = 10.8 µm. Concerning the potential dynamics, the shift in the peak time ∆ V

A-E is one order of magnitude smaller, going from 0.02 sec for L 2 = 0.9 µm to 0.03 sec for L 2 = 10.8 µm. The ratio r V

A-E of the electrical signal reaching (E) goes from 84% for L 2 = 0.9 µm, to 62% for L 2 = 10.8 µm.

Declarations

Authors' contributions: CG and SK designed the project. PP and SK built the code. CG and PP designed the application to neuroscience. PP run the simulations. CG, SK and PP wrote the manuscript.

Funding: Financial support from the EUR Spectrum and EUR Life in Université Côte d'Azur is gratefully acknowledged.

Competing interest statement: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.