Claire Alamichel

Juan Calvo

Erwan Hingant

Saoussen Latrach

Nathan Quiblier

Romain Yvinec

Signaling Pathway

Modeling compartmentalization within intracellular signaling pathway

come

Introduction

G Protein Coupled Receptors (GPCR) are a large class of transmembrane receptors: they are proteins located at the cell surface (plasma membrane) that are specialized in receiving (binding to) extracellular molecules, called ligand. The ligand-bound receptor will subsequently affect biochemical changes through the plasma membrane towards intracellular molecules. GPCR thus play a key role in the signal transduction, which allows long range cell communications mediated by ligand (e.g. hormones, cytokines, growth factors etc) and leads to major changes in the metabolism and activity of a cell [START_REF] Lefkowitz | A Brief History of G-Protein Coupled Receptors (Nobel Lecture)[END_REF][START_REF] Roth | G Protein-Coupled Receptor Signalling Networks from a Systems Perspective[END_REF]. In addition, GPCR form an important class of targeted pharmaceutical agents in many different physiological contexts. A receptor-activated signalling pathway, or biochemical cascade, is a chain of biochemical events activated by a receptor upon ligand binding. A series of signal transducers, or effector molecules, are activated within a reaction network (with typically many feedbacks and/or feedforward loops) and ultimately leads to a cell response, e.g. changes in gene expression, cell growth, metabolism [START_REF] Kholodenko | Cell-signalling dynamics in time and space[END_REF]. Recently, it has been shown that GPCR are pleiotropic: they are able to selectively activate different signalling pathways. Distinct biochemical events can indeed be modulated by the nature of the ligand (either native physiological ligand or pharmaceutical agent) and the specific 3D conformational structured adopted by the ligand-bound receptor. This pleiotropy sheds light into important cell regulation mechanisms that need to be taken into account to design efficient therapeutic strategies [START_REF] Kenakin | Biased Receptor Signaling in Drug Discovery[END_REF]. GPCR do not permanently stay at the plasma membrane. Rather, an important desensitization mechanism of signaling pathways is receptor endocytosis. Receptor endocytosis is the internalization inside the cell of a vesicle surrounded by an area of cell membrane which contains receptor (and other) molecules. Endocytosis is an active form of transport of molecules. The vesicles, called endosomes, and their molecular content, can further be degraded by cell machinery through hydrolytic enzymes, or be partly recycled back to the plasma membrane. These spatial movements induce, upon ligand binding, an heterogeneous population of receptors spread out between plasma membrane and a dynamic population of endosomes with distinct micro-environment. Within the context of receptor-activated signalling pathways, the processes of endocytosis, degradation and/or recycling are commonly referred to receptor trafficking, and provide further cell regulation mechanisms of signaling pathways [START_REF] Birtwistle | Endocytosis and signalling: A meeting with mathematics[END_REF]. Key evidence have indeed shown that the intra-cellular traffic of internalized receptors has a major impact in cell response to a given stimuli [START_REF] Jean-Alphonse | Regulation of GPCR signal networks via membrane trafficking[END_REF][START_REF] Vilardaga | Endosomal generation of cAMP in GPCR signaling[END_REF]. In this work, we will focus our applications on a specific class of receptors that share in common the activation of the same effector molecule, the cyclic adenosine monophosphate (cAMP). For a number of GPCR among this class, including the Beta-2 Adrenergic Receptors (β2AR) [START_REF] Kim | Spatiotemporal Characterization of GPCR Activity and Function during Endosomal Trafficking Pathway[END_REF], the parathyroid hormone receptor (PTHR) [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF], the luteinizing hormone receptor (LHR) [START_REF] Lyga | Persistent cAMP Signaling by Internalized LH Receptors in Ovarian Follicles[END_REF], or the follicle-stimulating hormone receptor (FSHR) [START_REF] Sayers | Intracellular Follicle-Stimulating Hormone Receptor Trafficking and Signaling[END_REF], after ligand binding, the production of cAMP occurs first at the plasma membrane and later on from a highly dynamic pool of endosomes following internalization of the receptor by endocytosis. This spatio-temporal dimension of signaling has been found to have a significant impact on physiological functions, such as the control of serum calcium by PTHR signaling [START_REF] White | Spatial bias in cAMP generation determines biological responses to PTH type 1 receptor activation[END_REF], or the resumption of meiosis by LHCGR signaling [START_REF] Lyga | Persistent cAMP Signaling by Internalized LH Receptors in Ovarian Follicles[END_REF]. Thus, to faithfully represent the complexity of signalling pathways, we need to take into account the dynamic of the transient pool of specialised endosomes [START_REF] Birtwistle | Endocytosis and signalling: A meeting with mathematics[END_REF][START_REF] Sorkin | Endocytosis and signalling: Intertwining molecular networks[END_REF][START_REF] Villaseñor | Signal processing by the endosomal system[END_REF] following receptor stimulation, and its role on the reaction networks involved in the signalling pathways. The current biological hypothesis is that the endosomal compartments provide a dynamic and heterogeneous compartmentalised structure that allows specialised effector molecules to be separated from the bulk cytoplasm (physically separated through a lipid bilayer) in order to have a proper function of the cell response. Classical ways to model the dynamic behaviour of signalling pathways use Chemical Reaction Networks [START_REF] Ingalls | Mathematical Modeling in Systems Biology[END_REF], either in a deterministic formalism using ordinary differential equations [START_REF] Feinberg | Foundations of Chemical Reaction Network Theory[END_REF], or in a stochastic formalism (typically when few molecules are present) using continuous-time Markov chains [START_REF] Anderson | Stochastic Analysis of Biochemical Systems[END_REF]. Both approaches typically assume the law of mass action and an idealised homogeneous environment. When spatial dynamics is important to take into account, one may use reaction-diffusion models to represent for instance spatial gradients [START_REF] Kholodenko | Cell-signalling dynamics in time and space[END_REF], or compartmental models, to physically represent segregation between static compartments [START_REF] Weddell | Integrative meta-modeling identifies endocytic vesicles, late endosome and the nucleus as the cellular compartments primarily directing RTK signaling[END_REF]. Note that in most applied literature, like in the field of epidemiology (e.g. SIR model), compartmental models refers to finite-dimensional dynamical systems where one variable is associated to each compartment, and the model described fluxes between compartments. Up to our knowledge, relatively few works have addressed the issue of representing explicitly the segregation of molecules in a dynamic environment. The peculiarity of receptor trafficking at play within signalling pathways is that the endosomal compartments are created upon receptor activation, and their number, their size and their molecular content evolved dynamically, within a similar time scale than that of the signalling pathways activation. A first attempt of defining such models, within the context of signalling pathways, dates back to [START_REF] Foret | A general theoretical framework to infer endosomal network dynamics from quantitative image analysis[END_REF] and uses a deterministic population dynamics formalism to follow a population of compartments, structured by their size and molecular content, and which undergo coagulation-fragmentation like dynamics, representing endocytosis, fusion, fission, recycling and degradation. Recently, a stochastic counterpart has been proposed by [START_REF] Duso | Stochastic reaction networks in dynamic compartment populations[END_REF][START_REF] Pietzsch | Compartor: a toolbox for the automatic generation of moment equations for dynamic compartment populations[END_REF][START_REF] Anderson | Stochastic Reaction Networks Within Interacting Compartments[END_REF]. To the best of our knowledge, a model that represents the chemical reactions that take place both at plasma membrane and within a dynamic population of endosomal compartment has never been considered and represents a novelty of our current work. In this work, we define and study the long-time behavior and numerical schemes of minimal deterministic models that can represent compartmentalised signalling pathways, taking inspiration from [START_REF] Foret | A general theoretical framework to infer endosomal network dynamics from quantitative image analysis[END_REF]. Furthermore, we provide two simple examples of our model that show that this model is able to represent qualitatively main experimental observations on cAMP signaling from [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF][START_REF] White | Spatial bias in cAMP generation determines biological responses to PTH type 1 receptor activation[END_REF]. In section 2, we describe two minimal models, structured with respectively one or two variables, that can represent the size-distribution of the endosomal compartment population and their molecular content. In section 3, we study the long time behavior of the model structured with a single variable. In section 4, we present a numerical scheme for the more general model, structured with two variables. In section 5, we present some numerical simulations that provide qualitative comparisons with experimental observations.

Modelling compartmentalised signalling pathways

The first objective is to be able to define a model to simulate the endosomal compartment dynamics from their size structure perspective only. We adopt a deterministic population dynamic approach, where individuals are structured by a single positive variable (their size). From biological observations, the main processes that shape the size distribution of the endosomal population include:

• Endocytosis: creation of a de-novo compartment from the cell membrane;

• Removal of compartment: either recycling back to the cell membrane, or degradation through lysosomal pathways; • Fusion: binary coagulation of compartments. Let g = g(t, r) be the population density of endosomal compartments at time t ∈ R + and size (volume) r ∈ R + . The evolution equation for g, that takes into account the three mechanisms above is given by, for all t > 0 and r > 0:

∂g(t, r) ∂t = Q(g, g)(t, r) coagulation + α(r) endocytosis -γ(r)g(t, r) removal , (M1)
where

Q(g, g)(t, r) = 1 2 r 0 κ(r -r ′ , r ′ )g(t, r -r ′ )g(t, r ′ )dr ′ - ∞ 0 κ(r, r ′ )g(t, r)g(t, r ′ )dr ′ . (1)
In the sequel, we denote g(0, r) = g 0 (r) the initial condition of Eq. (M1). A similar equation may be found in [START_REF] Foret | A general theoretical framework to infer endosomal network dynamics from quantitative image analysis[END_REF][START_REF] Alexandrov | Dynamics of intracellular clusters of nanoparticles[END_REF]. The coagulation operator is quadratic: two compartments of respective size r, r ′ fuse in one compartment of size r + r ′ at rate κ(r, r ′ ). One obtains Eq. ( 1) by considering all compartments that reach size r and all that leave size r. The kernel κ ≥ 0 is symmetric (κ(r, r ′ ) = κ(r ′ , r)), such that the coagulation operator Eq. ( 1) preserves mass ( rQ(g, g)(r)dr = 0). In particular, we have (at least formally),

d dt ∞ 0 rg(t, r)dr = ∞ 0 rα(r)dr - ∞ 0 γ(r)rg(t, r)dr . (2)
Endocytosis is a zero-order process at rate α (source term), compartment removal is a first-order process and occurs at a rate γ. In the sequel, model (M1) will be referred as our 1D model.

The second modelling step is to include molecular content into Eq. (M1). One may first think as an additional one dimensional structure variable 1 . The first example we have in mind is the quantity of (active) receptor within each compartment, which is of primary interest to represent receptor trafficking within cells (see application in subsection 5.1). The second example we have in mind is the production of cAMP effector molecules, both at the plasma membrane and within each endosomes (see application in subsection 5.2). Thanks to this second structuring variable, besides compartment dynamics, we aim to represent:

• Biochemical reactions inside each compartment. Reaction rates are dependent on local abundances of molecular species, as well as the size (and more generally other physical variable like pH) of the compartment. • Biochemical reactions that occur at the plasma membrane.

• Molecular conservation laws that may hold at each event that modifies the compartments, in particular between membrane and compartments. Hence we will now distinguish between compartment degradation and compartment recycling. Let f = f (t, r, a) be the population density of compartments at time t ∈ R + , size (volume) r ∈ R + and molecular content a ∈ R + . Let also M = M (t) ∈ R + be the molecular quantity at the plasma membrane. The joint evolution equation for f, M is , for all t > 0, r > 0 and a > 0,

∂f (t, r, a) ∂t + ∂ (V (r, a)f (t, r, a)) ∂a reactions = Q(f, f )(t, r, a) coagulation + α(r, a, M (t)) endocytosis -γ(r, a)f (t, r, a) degradation -λ(r, a)f (t, r, a)
recycling , (M2a) 1 In future work, this additional structure variable could be an arbitrary finite dimensional variable to represent other molecular actors of signalling pathways that are either physically located at the plasma membrane, in the endosomal compartments or at the vicinity of those.

dM (t) dt = J M (M (t)) reactions - ∞ 0 ∞ 0 aα(r, a, M (t))dadr endocytosis + ∞ 0 ∞ 0 aλ(r, a)f (t, r, a)dadr recycling , (M2b)
where

Q(f, f )(t, r, a) = 1 2 r 0 a 0 κ(r -r ′ , r ′ )f (t, r -r ′ , a -a ′ )f (t, r ′ , a ′ )da ′ dr ′ - ∞ 0 ∞ 0 κ(r, r ′ )f (t, r, a)f (t, r ′ , a ′ )da ′ dr ′ . ( 3
)
In the sequel, we denote f (0, r, a) = f 0 (r, a) and M (0) = M 0 the initial conditions of Eqs. (M2a)-(M2b).

In Eq. (M2a), the extension of the coagulation and endocytosis processes to include the second structuring variable a is clear. Note that we chose the coagulation kernel to be dependent on the size of the compartments (not their molecular content), for the sake of simplicity. The removal terms are now splitted in two: the degradation occurs at rate γ, and the recycling occur at rate λ. The transport term represents the biochemical reactions that modify the molecular content within each compartment, and which occur at rate V (r, a). Also, the endocytosis rate α is necessarily dependent on the molecular content of the plasma membrane M , to avoid negative values for M , e.g. α(r, a, 0) = 0. In Eq. (M2b), J M represent biochemical reactions that occur at the plasma membrane, and the two integral terms represent the molecules that are lost or gained at the plasma membrane through respectively the endocytosis or recycling processes. The Eq. (M2b) on the scalar variable M was not present in [START_REF] Foret | A general theoretical framework to infer endosomal network dynamics from quantitative image analysis[END_REF] and, to the best of our knowledge, it is a novelty of our model. Still interpreting the molecular content as a quantity of receptors, it allows to represent conservation laws of receptors between plasma membrane and endosomal compartments, giving a satisfactory representation of the biological concept of receptor trafficking: receptors (together with other molecules) undergo directed movement back-and-forth between plasma membrane and intra-cellular endosomal compartments. In particular, we have (at least formally),

d dt ∞ 0 ∞ 0 rf (t, r, a)drda = ∞ 0 ∞ 0 rα(r, a, M (t))drda - ∞ 0 ∞ 0 γ(r, a)rf (t, r, a)drda - ∞ 0 ∞ 0 λ(r, a)rf (t, r, a)drda , (4)
and

d dt ∞ 0 ∞ 0 af (t, r, a)drda + M (t) = J M (M (t)) + ∞ 0 ∞ 0 V (r, a)f (t, r, a)drda - ∞ 0 ∞ 0 γ(r, a)af (t, r)drda . ( 5
)
The choice of the rate functions α, γ, λ and κ mainly governs the compartment size dynamics. In [START_REF] Foret | A general theoretical framework to infer endosomal network dynamics from quantitative image analysis[END_REF][START_REF] Alexandrov | Dynamics of intracellular clusters of nanoparticles[END_REF], the authors consider a one-dimensional model (M1) with constant coagulation kernel, constant first-order rate γ, and exponentially decaying source rate α(r) = Ce -r/a0 . The choice of the rate functions V and J M will mainly impact the molecular content at the plasma membrane and within compartments. Its choice may be guided by the underlying chemical reactions that take place at the plasma membrane and within compartments. Linear or polynomial functions may thus be suitable to represent mass action kinetic law.

In the sequel, the model given by Eqs. (M2a)-(M2b) will be referred as our 2D model.

In the remaining, in Eq. (M1), Eqs. (M2a)-(M2b), all rate functions and initial conditions are nonnegative. Further, we do not address in this paper the well-posedness of solutions of either Eq. (M1) or Eqs. (M2a)-(M2b) and rather take for granted that a unique sufficiently smooth nonnegative solution is given.

Long time behavior of the 1D model (M1)

In this section we provide sufficient conditions so that Eq. (M1) exhibits a unique globally stable steady state.

The large-time behaviour of this equation -also known as coagulation equation with source and efflux [START_REF] Chae | Existence and uniqueness for spatially inhomogeneous coagulation equation with sources and effluxes[END_REF], has been studied first in [START_REF] Gajewski | On a first order partial differential equation with nonlocal nonlinearity[END_REF] with drift, in [START_REF] Dubovskii | Mathematical Theory of Coagulation[END_REF][START_REF] Laurençot | Stationary solutions to smoluchowski's coagulation equation with source[END_REF] without efflux, in [START_REF] Vásquez | Ecuaciones en derivadas parciales para el análisis de modelos biopoliméricos[END_REF] with a bounded coagulation kernel and in [START_REF] Ghosh | On equilibrium solution to a singular coagulation equation with source and efflux[END_REF] with a singular coagulation kernel. Here, we limit ourselves to give a self-contained proof of exponential stability of the steady state in L 1 with bounded coagulation kernel. We mainly use a contraction argument, taking inspiration from [START_REF] Collet | Lifshitz-slyozov equations: The model with encounters[END_REF].

Theorem 1. Let g 0 , α, κ, γ nonnegative. Assume α is integrable, κ and γ are bounded and morevover

inf γ = γ 0 > 0. If 3∥κ∥ L ∞ ∥α∥ L 1 < γ 2 0 , (6)
then there exists a unique nonnegative stationary solution in L 1 (R + ) of Eq. (M1), denoted by g ∞ . Moreover,

∥g ∞ ∥ L 1 ≤ ∥α∥ L 1 γ0
, and for every solution

g ∈ C(R + , L 1 (R + )) we have lim t→+∞ ∥g(t) -g ∞ ∥ L 1 = 0 .
The convergence is at least exponential with rate

γ 2 0 -3∥κ∥ L ∞ ∥α∥ L 1 γ0 > 0.
We do not expect condition (6) to be optimal. See the discussion in Section 4.2. We assume here that wellposedness of a unique nonegative solution g ∈ C(R + , L 1 (R + )) of Eq. (M1) is given. We note that this implicitly imposes conditions on the rate functions and the initial condition, which will at least need to be nonnegative and integrable. See [START_REF] Canizo | Some problems related to the study of interaction kernels: coagulation, fragmentation and diffusion in kinetic and quantum equations[END_REF] for more details.

Proof. First we prove existence and uniqueness of a stationary solution g ∞ thanks to a Banach fixed point argument. Let us define

X = g ∈ L 1 (R + ) : g ≥ 0, ∥g∥ L 1 ≤ ∥α∥ L 1 γ 0 . Consider a constant K ≥ ∥γ∥ L ∞ + ∥κ∥ L ∞ ∥α∥ L 1 γ0 .
For g ∈ X we define

T K g = 1 K (α -γg + Q(g, g) + Kg) .
The operator T K is well-defined since Q(g, g) is also well-defined for any integrable g provided that κ is bounded [30, Lemma 3]. We aim to apply Banach's fixed point theorem for T K on X -which is a closed subset of the Banach space L 1 (R + ); this will give a stationary solution of Eq. (M1). Let g ∈ X, we have that

T K g ≥ 1 K (K -∥γ∥ L ∞ -∥κ∥ L ∞ ∥g∥ L 1 ) g . (7)
Indeed, g and α are positive and κ is bounded, thus

Q(g, g) ≥ -g(r) ∞ 0 κ(r, r ′ )g(r ′ )dr ′ ≥ -∥κ∥ L ∞ ∥g∥ L 1 g(r) ,
which leads to Eq. ( 7) thanks to the boundedness of γ. The fact that ∥g∥

L 1 ≤ ∥α∥ L 1 γ0
and the condition on K entail the positivity of T K f . Then, because

∞ 0 Q(g, g) dr = - 1 2 ∞ 0 ∞ 0 κ(r ′ , r)g(r ′ )g(r)dr ′ dr ≤ 0 ,
we deduce that

∥T K g∥ L 1 = ∞ 0 1 K (α(r) -γg(r) + Q(g, g)(r) + Kg(r)) dr ≤ ∥α∥ L 1 K + K -γ 0 K ∥g∥ . But K ≥ ∥γ∥ L ∞ ≥ γ 0 , thus, for g ∈ X, K∥T K g∥ L 1 ≤ ∥α∥ L 1 + (K -γ 0 ) ∥α∥ L 1 γ 0 ≤ K ∥α∥ L 1 γ 0 ,
and we conclude that T K g belongs to X. We now show that T K is a contraction on X. Let g and h in X. It is straightforward to check that

∥T K g -T K h∥ L 1 ≤ (1 - γ 0 K )∥g -h∥ L 1 + 1 K ∥Q(g, g) -Q(h, h)∥ L 1 .
By a simple computation (see e.g. [START_REF] Collet | Lifshitz-slyozov equations: The model with encounters[END_REF]Lemma 3] or [START_REF] Vásquez | Ecuaciones en derivadas parciales para el análisis de modelos biopoliméricos[END_REF][START_REF] Canizo | Some problems related to the study of interaction kernels: coagulation, fragmentation and diffusion in kinetic and quantum equations[END_REF]),

∥Q(g, g) -Q(h, h)∥ L 1 ≤ 3 2 ∥κ∥ L ∞ (∥g∥ L 1 + ∥h∥ L 1 )∥g -h∥ L 1 ≤ 3∥κ∥ L ∞ ∥α∥ L 1 γ 0 ∥g -h∥ L 1 . (8)
Thus, we have

∥T K g -T K h∥ L 1 ≤ 1 + 1 K (3∥κ∥ L ∞ ∥α∥ L 1 γ 0 -γ 0 ) ∥g -h∥ L 1 .
Then the hypothesis Eq. ( 6) allows us to conclude that there exists a unique fixed point to T K in X. Moreover, this is the unique stationary solution which is positive and belongs to L 1 . Indeed, assume we have a positive stationary solution g /

∈ X i.e. in L 1 satisfying γ 0 g > ∥α∥ L 1 . Then 0 = ∞ 0 (α(r) -γ(r)g(r) + Q(g, g)(r)) dr ≤ ∥α∥ L 1 -γ 0 ∞ 0 g(r) dr < 0 ,
which is a contradiction. We now turn to the proof of asymptotic stability. Let g ∈ C(R + , L 1 (R + )) be a nonnegative solution to Eq. (M1) in the sense of distributions and let g ∞ ∈ L 1 (R + ) be the stationary solution. We first provide a bound on g. We have

d dt ∞ 0 g(t, r)dr = ∥α∥ L 1 -γ 0 ∞ 0 g(t, r)dr + ∞ 0 Q(g, g)(t, r)dr ≤ ∥α∥ L 1 -γ 0 ∞ 0 g(t, r)dr .
Thus,

∥g(t, •)∥ L 1 ≤ ∥g(0, •)∥ L 1 e -γ0t + ∥α∥ L 1 γ 0 .
It is a classical computation that

∂ ∂ t |g -g ∞ | = -γ(r)|g -g ∞ | + (Q(g, g) -Q(g ∞ , g ∞ )) sign(g -g ∞ ).
Thus, using Eq. ( 8),

d dt ∥g(t, •) -g ∞ ∥ L 1 ≤ -γ 0 + 3 2 ∥κ∥ L ∞ ∥g(0, •)∥ L 1 e -γ0t + ∥α∥ L 1 γ 0 + ∥g ∞ ∥ L 1 ∥g(t, •) -g ∞ ∥ L 1 ≤ 3∥κ∥ L ∞ ∥α∥ L 1 γ 0 -γ 0 + 3 2 ∥κ∥ L ∞ ∥g(0, •)∥ L 1 e -γ0t ∥g(t, •) -g ∞ ∥ L 1 .
We conclude that

∥g(t, •) -g ∞ ∥ L 1 ≤ ∥g(0, •) -g ∞ ∥ L 1 e 3∥κ∥ L ∞ 2γ0 ∥g(0,•)∥ L 1 - γ 2 0 -3∥κ∥ L ∞ ∥α∥ L 1 γ0 t ,
which ends the proof. □

Numerical scheme for 1D and 2D models

In this section, we detail our Finite Volume numerical scheme used to simulate our models, and numerically illustrate their properties. In order to tackle the complexity of the model, we decided to construct a conservative scheme that preserves moments of the solution. To that, we use of the finite volume scheme framework, after reformulating the coagulation operator in a divergence form in the spirit of [START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF][START_REF] Hingant | Derivation and mathematical study of a sorption-coagulation equation[END_REF]. The transport term is treated by an up-wind scheme. Remaining terms are approximated in a standard way. Time is treated by a first order Euler explicit approximation.

Finite volume scheme

We want to write a numerical scheme of the 2D model, Eqs. (M2a)-(M2b) using a finite volume method.

First we can remark that, for r > 0 and a > 0, we can rewrite our Eq. (M2a) in the following conservative form [START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF] ∂

t f (t, r, a) = 1 ra ∂ r ∂ a C (f ) (t, r, a) -∂ a (V (r, a) f (t, r, a)) + α (r, a, M (t)) -λ (r, a) f (t, r, a) -γ (r, a) f (t, r, a) , (9)
where

C (f ) (t, r, a) = r 0 a 0 r ′ a ′ Q (f, f ) (t, r ′ , a ′ )da ′ dr ′ .
From Eq. ( 9), we will first detail the truncation we use (step (i)) and the discretization (step (ii)). Then, we detail the finite volume approximation of the coagulation and first-order transport operator in the right-hand side of Eq. ( 9) (step (iii)). We finally sum-up the numerical scheme in the last step (step (iv)).

(i) Truncation. As in [START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF][START_REF] Hingant | Derivation and mathematical study of a sorption-coagulation equation[END_REF], to study our equation, we will truncate the size variable to a maximal value R > 0 and the quantity of reactants variable to a maximal value A > 0 and we will choose a truncation of the functional C. We chose the following truncation, given by, for any r ≤ R and a ≤ A,

C RA c (f ) (t, r, a) = 1 2 r 0 a 0 r ′ a ′ r ′ 0 a ′ 0 κ (r ′ -r ′′ , r ′′ ) f (t, r ′ -r ′′ , a ′ -a ′′ ) f (t, r ′′ , a ′′ ) da ′′ dr ′′ da ′ dr ′ - r 0 a 0 r ′ a ′ f (t, r ′ , a ′ ) R-r ′ 0 A-a ′ 0 κ (r ′ , r ′′ ) f (t, r ′′ , a ′′ ) da ′′ dr ′′ da ′ dr ′ . (10)
The truncation in Eq. ( 10) ensures that no cluster of size larger than R and molecular content larger than A arise. It can be obtained from Eq. (3) using the truncated kernel κ(r, r ′ )1 r+r ′ <R 1 a+a ′ <A instead of κ(r, r ′ ). The equation we will numerically approximate is a truncated version of Eq. ( 9) on the time interval [0, T ], where T > 0. We then look at the equation, for any 0 < t < T , 0 < r ≤ R and 0 < a ≤ A,

∂ t f RA (t, r, a) = 1 ra ∂ r ∂ a C RA c f RA (t, r, a) -∂ a V (r, a) f RA (t, r, a) + α (r, a, M (t)) -λ (r, a) f RA (t, r, a) -γ (r, a) f RA (t, r, a) . (11)
Heuristically, we expect f RA to be close to f as long as the mass of f outside [0, R/2] × [0, A/2] is small (so that coagulation that leads to compartment of size larger than R and molecular content larger than A are unlikely). We call Eq. ( 11) a conservative truncation of Eq. ( 9) because the choice of the truncation of the coagulation operator does not lead to a loss of mass in the first moments. In particular, f RA satisfies, at least formally, (compare to Eqs. ( 4)-( 5)),

d dt ∞ 0 ∞ 0 rf RA (t, r, a)drda = ∞ 0 ∞ 0 rα(r, a, M )drda - ∞ 0 ∞ 0 γ(r, a)rf RA (t, r, a)drda - ∞ 0 ∞ 0 λ(r, a)rf RA (t, r, a)drda , (12)
and

d dt ∞ 0 ∞ 0 af RA (t, r, a)drda + M (t) = J M (M (t)) + ∞ 0 ∞ 0 V (r, a)f RA (t, r, a)drda - ∞ 0 ∞ 0 γ(r, a)af RA (t, r)drda . ( 13
)
We choose to use a conservative truncation for our scheme in order to construct a scheme that preserves the conservation properties ( 12)-( 13), as we will verify it in the case of pure coagulation.

(ii) Grids definition. Let I r ∈ N. We discretize the size interval [0, R] into I r intervals. We denote by

r i-1 2 i∈{1,...,I r +1}
a regular mesh of [0, R] with size step ∆r and we set

r i = r i-1 2 + r i+ 1 2 2 = i - 1 2 ∆r, i ∈ {1, . . . , I r } .
Let I a ∈ N. We discretize the size interval [0, A] into I a intervals. We denote by a j-1 2 j∈{1,...,I a +1} a regular mesh of [0, A] with step ∆a and we set

a j = a j-1 2 + a j+ 1 2 2 = j - 1 2 ∆a, j ∈ {1, . . . , I a } .
For all i ∈ {1, . . . , I r } and j ∈ {1, . . . , I a }, we set

Λ ij = r i-1 2 , r i+ 1 2 × a j-1 2 , a j+ 1 2 .
Let ∆t > 0 be the time step. We discretize [0, T ] by the set of points {t n = n ∆t, n ∈ {0, . . . , N }}, where N = T ∆t .

(iii) Finite volume approximation. For all n ∈ {0, . . . , N }, i ∈ {1, . . . , I r } and j ∈ {1, . . . , I a }, we denote f n i,j an approximation of the function f RA at the point (t n , r i , a j ). We will recursively calculate f n i,j such that

f n i,j ≈ 1 ∆r∆a Λij f RA (t n , r, a) dadr.
Similarly, we denote by M n an approximation of M (t n ). First, we write the explicit forward Euler scheme in time associated to the Eq. ( 11). For all n ∈ {0, . . . , N -1}, we have

f RA t n+1 , r, a -f RA (t n , r, a) ∆t = 1 ra ∂ r ∂ a C RA c f RA (t n , r, a) -∂ a V (r, a) f RA (t n , r, a) + α (r, a, M n ) -λ (r, a) f RA (t n , r, a) -γ (r, a) f RA (t n , r, a) +O(∆t). ( 14
)
For each i ∈ {1, . . . , I r } and j ∈ {1, . . . , I a }, we will integrate Eq. ( 14) over Λ ij . We first deal with the coagulation operator. Approximating ra by the constant r i-1 2 a j-1 2 on Λ ij , and integrating, we obtain

Λij 1 ra ∂ r ∂ a C RA c f RA (t n , r, a) dadr ≈ 1 r i-1 2 a j-1 2 (C RA c f RA t n , r i+ 1 2 , a j+ 1 2 -C RA c f RA t n , r i-1 2 , a j+ 1 2 -C RA c f RA t n , r i+ 1 2 , a j-1 2 + C RA c f RA t n , r i-1 2 , a j- 1 2
). ( 15)

Using a change of variable r ′ -r ′′ → r ′ in Eq. ( 10), and decomposing the first two integrals in telescopic sums we have,

C RA c f RA t n , r i+ 1 2 , a j+ 1 2 = 1 2 i k=1 j m=1 r k+ 1 2 r k-1 2 a m+ 1 2 a m-1 2 r i+ 1 2 -r ′ 0 a j+ 1 2 -a ′ 0 a ′ + a ′′ r ′ + r ′′ κ r ′ , r ′′ f RA t n , r ′ , a ′ f RA t n , r ′′ , a ′′ da ′′ dr ′′ da ′ dr ′ - i k=1 j m=1 r k+ 1 2 r k-1 2 a m+ 1 2 a m-1 2 r ′ a ′ f RA t n , r ′ , a ′ R-r ′ 0 A-a ′ 0 κ r ′ , r ′′ f RA t n , r ′′ , a ′′ da ′′ dr ′′ da ′ dr ′ . Further, setting κ k,k ′ = 1 ∆r 2 r k+ 1 2 r k-1 2 r k ′ + 1 2 r k ′ -1 2
κ (r, r ′ ) dr ′ dr, and approximating again product terms ra by their left value on Λ ij , we approximate C RA c by:

C RA c f RA t n , r i+ 1 2 , a j+ 1 2 ≈ 1 2 (∆r∆a) 2 i k=1 j m=1 i-k+1 k ′ =1 j-m+1 m ′ =1 a m-1 2 + a m ′ -1 2 r k-1 2 + r k ′ -1 2 κ k,k ′ f n k,m f n k ′ ,m ′ -(∆r∆a) 2 i k=1 j m=1 I r -k+1 k ′ =1 I a -m+1 m ′ =1 a m-1 2 r k-1 2 κ k,k ′ f n k,m f n k ′ ,m ′ , (16)
Plugging this latter approximation (16) into Eq. ( 15), we finally define the approximation of the coagulation operator, as

Q n i,j (f, f ) = (∆r∆a) 2 1 2 i k=1 j m=1 κ k,i-k+1 f n k,m f n i-k+1,j-m+1 -f n i,j I r -i+1 k=1 I a -j+1 m=1 κ i,k f n k,m . (17)
Note that the approximation Q n i,j of Q f RA , f RA is conservative (in the sense that it preserves the first moments in r and a), as we will show in the next section. Then, we use an upwind approximation for the transport term, together with an integration of transport velocity over r using the midpoint rule:

Λij ∂ a V (r, a) f RA (t, r, a) dadr = r i+ 1 2 r i-1 2 V r, a j+ 1 2 f RA t, r, a j+ 1 2 dr - r i+ 1 2 r i-1 2 V r, a j-1 2 f RA t, r, a j-1 2 dr ≈ W i,j+ 1 2 f RA t, r i , a j+ 1 2 ∆r -W i,j- 1 2 f RA t, r i , a j-1 2 ∆r ≈ ∆r • A up W i,j+ 1 2 , f n i,j , f n i,j+1 -A up W i,j- 1 2 , f n i,j-1 , f n i,j , (18)
where

W i,j-1 2 = 1 2 V r i+ 1 2 , a j-1 2 + V r i-1 2 , a j-1 2
, f n i,0 = f n i,I r +1 = 0, and the operator A up is

A up (u, f + , f -) = uf + if u ≥ 0, uf -if u < 0.
(iv) Finite volume scheme. Finally, the scheme of the model (M2a)-( M2b) is given by: • Initialization for i ∈ {1, . . . , I r } and j ∈ {1, . . . , I a } we set

f 0 i,j = 1 ∆r∆a Λij f (0, r, a) dadr.
• Time iteration: for all n ∈ {1, . . . , N }, i ∈ {1, . . . , I r } and j ∈ {1, . . . , I a }, using Eqs. ( 17)-( 18), and using an explicit forward Euler scheme in time for M , we set:

f n+1 i,j =f n i,j + ∆t ∆r∆a Q n i,j (f, f ) - ∆t ∆a A up W i,j+ 1 2 , f n i,j , f n i,j+1 -A up W i,j- 1 2 , f n i,j-1 , f n i,j + ∆t ∆r∆a Λij α (r, a, M n ) dadr -f n i,j Λij λ (r, a) dadr -f n i,j Λij γ (r, a) dadr . M n+1 =M n + ∆t J M (M n ) -∆t R 0 A 0 a ′ α (r ′ , a ′ , M n ) da ′ dr ′ + ∆t I r i=1 I a j=1 f n i,j Λij aλ (r, a) dadr. ( 19
)
The remaining integrals in Eq. ( 19) are calculated using an automatic adaptive numerical integration method, detailed in [START_REF] Genz | Remarks on algorithm 006: An adaptive algorithm for numerical integration over an Ndimensional rectangular region[END_REF] and implemented in the Julia package HCubature2 .

Remark 1. With the same tools, we can write a scheme in dimension 1 for model (M1), noticing that if we set h (t, r) = rg (t, r) we have

∂ t h (t, r) = -∂ r J (g) (t, r) + rα (r, M (t)) -γ (r) h (t, r) ,
where J (g) is defined as follows

J (g) (t, r) = r 0 ∞ r-r ′ r ′ κ (r ′ , r ′′ ) g (t, r ′ ) g (t, r ′′ ) dr ′′ dr ′ .
A conservative truncation is then given by the following:

J R c (g) (t, r) = r 0 R-r ′ r-r ′ r ′ κ (r ′ , r ′′ ) g (t, r ′ ) g (t, r ′′ ) dr ′′ dr ′ .

Numerical tests

In this section, we will investigate the behavior of the numerical scheme given by Eq. ( 19) in some particular cases. We aim to verify the conservation laws satisfied by the numerical scheme, the consistency and convergence properties of the numerical scheme, as well as to illustrate the long-time behavior we proved in Theorem 1.

Conservation laws and consistency. First, we will evaluate the numerical scheme in the case of a pure coagulation model (α = γ = λ = V = J M = 0). We consider the pure coagulation equation:

∂ t f (t, r, a) = Q (f (t) , f (t)) (r, a) , with t > 0, r > 0 and a > 0. ( 20
)
The discrete equation associated to Eq. ( 20) is given by, as a special case of Eq. ( 19)

when α = γ = λ = V = J M = 0, f n+1 i,j = f n i,j + ∆t ∆r∆a Q n i,j (f, f ) , n ∈ {1, . . . , N } , i ∈ {1, . . . , I r } and j ∈ {1, . . . , I a } , (21)
and where Q n i,j is defined in Eq. ( 17). We remark that we keep the same properties on the moments dynamics for the two Eqs. ( 20) and [START_REF] Pietzsch | Compartor: a toolbox for the automatic generation of moment equations for dynamic compartment populations[END_REF]. For some test function φ, we define H (φ, t) the moment of f associated with the function φ at time t > 0,

H (φ, t) := ∞ 0 ∞ 0 φ (r, a) f (t, r, a) dadr ,
and H n (φ) its discrete analogue moment, associated with the function φ at time t n with n ∈ {1, . . . , N },

H n (φ) := ∆r∆a I r i=1 I a j=1 φ r i-1 2 , a j-1 2 f n i,j . (22)
Then

d dt H (φ, t) = d dt ∞ 0 ∞ 0 φ (r, a) f (t, r, a) dadr = ∞ 0 ∞ 0 φ (r, a) Q (f (t) , f (t)) (r, a) dadr = 1 2 ∞ 0 ∞ 0 ∞ 0 ∞ 0 [φ (r + r ′ , a + a ′ ) -φ (r, a) -φ (r ′ , a ′ )] κ (r, r ′ ) f (t, r, a) f (t, r ′ , a ′ ) da ′ dr ′ dadr ,
and similarly

H n+1 (φ) -H n (φ) ∆t = (∆r∆a) 2 2 I r i=1 I a j=1 I r -i+1 k=1 I a -j+1 m=1 φ r i+k-1-1 2 , a j+m-1-1 2 -φ r i-1 2 , a j-1 2 -φ r k-1 2 , a m-1 2 κ i,k f n i,j f n k,m .
Thus choosing φ (r, a) = r or φ (r, a) = a, we have that the first-order moments are constant in time both at the discrete and continuous levels, consistently with Eqs. ( 4)-( 5) and Eqs. ( 12)-( 13) for the pure coagulation case. We also have that the zeroth-order moments are non-increasing functions of time in both cases.

For the pure coagulation model given by Eq. ( 20), choosing an affine kernel κ (r, r ′ ) = K 0 + K 1 (r + r ′ ) leads to a closed moment equation in the form of an ODE system. Indeed, in such case, it is easy to see that we have the following ODE system for the moments of order 0 and 1:

                   d dt H (1, t) = - 1 2 K 0 (H (1, t)) 2 -K 1 H (1, t) H (r, t) , d dt H (r, t) = 0, d dt H (a, t) = 0. (23)
Eq. ( 23) is of the form of a Bernoulli differential equation, and its analytical solution can be computed. The pure coagulation case with an affine kernel is thus an appropriate setting to compare analytical solutions with moments associated to our numerical scheme (21) calculated by Eq. ( 22). For this purpose, we introduce the following notation, for x > 0:

N x (µ, σ) = 1 √ 2πσ e -1 2 ( x-µ σ ) 2 1 x>0 . (24)
We approximate the pure coagulation Eq. ( 20) using the scheme Eq. ( 21) with the following numerical parameters:

T = 1 , ∆t = 10 -4 , R = A = 10 , ∆r = ∆a = 0.25 . (25)
We study two cases. The first one is the case of a constant kernel with the following model parameters:

f 0 (
The second one is the case of an affine kernel with the following model parameters:

f 0 (r, a) = 0.5 • [N r (1.5, 0.15) × N a (0.5, 0.3) + N r (0.5, 0.3) × N a (1.5, 0.15)] , κ (r, r ′ ) = 0.5 + 0.1 (r + r ′ ) . (27)
As expected, we recover from our numerical scheme (21) that the moment of order 0 is a nonincreasing function and the moments of order 1 are constant (Figure 1). Moreover, we also observed that these numerical results are very close to the analytical solutions directly computed from Eq. ( 23) (Figures 2 and3). Relative error of moments of order 0 are increasing through time, as expected from the fact that the size-truncation of the numerical scheme ( 21) implies more and more error as compartments gets bigger in pure-coagulation dynamics.

Convergence. We now investigate the convergence property of our numerical scheme given by Eq. ( 19) in a more general setting. We want to illustrate the convergence of the numerical scheme as the size steps ∆r and ∆a decreases to 0. Regarding stability, we expect the numerical scheme to be conditionally stable for small enough ∆t ∆r∆a as our scheme is a first-order finite volume scheme, but the proof of stability is out of the scope of this paper. See for instance [START_REF] Bourgade | Convergence of a finite volume scheme for coagulation-fragmentation equations[END_REF] for stability results on the coagulation part. We thus fix a small enough time step ∆t and choose ∆r = ∆a = h with linearly decreasing h in a log2 scale, namely h m = h0 2 m for integers m from 0 to 11 and h 0 = 1 2 . The finite volume size is then decreasing from h 2 0 to h0

2 11
2 . We choose the finest solution as reference solution and compare the solutions defined on the coarser grids to this reference solution with the following different norms, defined for a function f and k, l ∈ N as follows:

∥f ∥ k,l = R 0 A 0 r k a l |f (r, a)| dadr.
Figure 1. Time evolution of the moments of order 0 and 1 (see inserted legend) computed from our pure coagulation numerical scheme [START_REF] Pietzsch | Compartor: a toolbox for the automatic generation of moment equations for dynamic compartment populations[END_REF], with the constant kernel and the model parameters [START_REF] Dubovskii | Mathematical Theory of Coagulation[END_REF] for the left picture and with the constant kernel and the model parameters [START_REF] Laurençot | Stationary solutions to smoluchowski's coagulation equation with source[END_REF] for the right one. The results are obtained with the numerical parameters [START_REF] Gajewski | On a first order partial differential equation with nonlocal nonlinearity[END_REF]. In both cases, the moment of order 0 is a nonincreasing function and the moments of order 1 are constant. With the choice of our initial condition, both moments of order 1 in r and a superimposed.

2. error of the moments of order 0 (left panel) and 1 (right panel) between moments computed from our pure coagulation numerical scheme [START_REF] Pietzsch | Compartor: a toolbox for the automatic generation of moment equations for dynamic compartment populations[END_REF] and from the analytical solution of the ODE system Eq. ( 23), in the case of a constant kernel [START_REF] Dubovskii | Mathematical Theory of Coagulation[END_REF] with the same parameters as in the left panel of Figure 1.

Those moment-based norms are natural within the context of coagulation equations, see e.g. Eqs. ( 12)- [START_REF] Sorkin | Endocytosis and signalling: Intertwining molecular networks[END_REF].

The fact that the grids are nested allows to compute the following errors, for m ∈ {0, ...10},

∥E(h m , t n )∥ k,l = h 2 11 R/h11 i=1 A/h11 j=1 r k i-1 2 a l j-1 2 | f n,hm i,j -f n,h11 i,j | , (28)
where the value of f n,hm i,j are extended from the coarser grid to the finer grid by taking its value constant over the subdomains defined by the finer grid. We compute the errors defined in Eq. ( 28) up to second order moments Figure 3. Relative error of the moments of order 0 (left panel) and 1 (right panel) between moments computed from our pure coagulation numerical scheme [START_REF] Pietzsch | Compartor: a toolbox for the automatic generation of moment equations for dynamic compartment populations[END_REF] and from the analytical solution of the ODE system Eq. ( 23), in the case of an affine kernel [START_REF] Laurençot | Stationary solutions to smoluchowski's coagulation equation with source[END_REF] with the same parameters as in the right panel of Figure 1.

(0 ≤ k + l ≤ 2) in two different cases as we detail now. For that aim, we introduce the notation, for x > 0,

P x (x, ϵ) = x + ϵ -x x + ϵ 1 3 x x 2 3 . (29)
(i) Case 1. We approximate the pure coagulation case Eq. ( 20) (Figure 4) with the numerical parameters ∆t = 0.0005 , R = A = 3 , ∆r = ∆a = h m ,

and the model parameters

(ii) Case 2. We approximate the general case Eqs. (M2a)-(M2b) (Figure 5), with the model parameters

and the same numerical parameters as in Eq. [START_REF] Collet | Lifshitz-slyozov equations: The model with encounters[END_REF]. In both cases, we observe that the error decreases linearly (in log-log scale) when the size of the discretization step h does. For both cases, the order of the scheme appears to be 1. 31)) . We plot the error as a function of the grid size h, in log-log scale, using the six different norms ∥E(h m , t n )∥ k,l defined in Eq. (28) for 0 ≤ k + l ≤ 2 (see inserted legend).

We show in Figure 6 two cases. In the first one (left hand-side), we use γ (r) = √ 3.1 together with Eq. ( 34), for which the condition Eq. ( 6) of Theorem 1 is satisfied. In the second one (right hand-side), we use γ (r) = 0.7 together with Eq. ( 34), for which the condition Eq. ( 6) is not satisfied. In both cases, the curves for t = 10 and t = 50 are superimposed, and we observe that the solution seems to converge to a stationary state. Numerically, while we verify the conclusion of Theorem 1 holds true, it seems that the condition Eq. ( 6) is too restrictive since the scheme stays stable in a wider range of parameters such that the condition Eq. ( 6) is not satisfied. 33)- [START_REF] Genz | Remarks on algorithm 006: An adaptive algorithm for numerical integration over an Ndimensional rectangular region[END_REF]. On the left hand-side the condition Eq. ( 6) is satisfied (γ (r) = √ 3.1) whereas on the right hand-side the condition Eq. ( 6) is not satisfied (γ (r) = 0.7). In both cases, the curves for t = 10 and t = 50 are superimposed.

Applications

In this section we provide two applications of our 2D model given by Eqs. (M2a)-(M2b), that show good qualitative agreement with published experimental data on receptor-activated signaling pathways. Our simulations also provide additional insight, which calls for new experiments.

Receptor trafficking

The 2D model given by Eqs. (M2a)-( M2b) is ideally suited to model receptor trafficking from plasma membrane to endocytic comparments. As detailed in the introduction, the GPCR are typically located at the surface of the cells, on the plasma membrane. Upon ligand binding, the GPCR activate several signalling pathways as well as their own internalization through endocytosis. The vesicles still carry the internalised receptors on their surface. Endosomes are then sorted thanks to complex processes which are not yet fully understood, but that depends both on the nature of the receptor and the ligand. Internalized receptors can indeed commit to several endosomal compartments of different kinds, and be recycled at the cell surface, which could impact on the kinetic profile of the receptor and its signalling pathways. Consequently, endocytosis regulate receptor cell surface density and signaling profile, and endosomal targeting of receptors may produce specificity in the signaling pathways. In particular, it has been shown [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF] that the LHR and the B2AR are two GPCR that undergo divergent trafficking to distinct endosomal compartments. B2AR traffics mostly to early endosomes (EEs) and LHR to pre-early endosomes (pre-EEs). In short, the authors in [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF] demonstrate that LHR endosome sizes increased over time quickly before reaching a plateau, producing a small endosome population (400-500 nm of diameter). The mean B2AR endosome sizes are bigger (1200-1400 nm of diameter). They also see that B2AR was more internalized than LHR in percentage but both receptors are equally recycled. The objective of T = 30, ∆t = 10 -1 , I r = I a = 30, R = A = 2000.

(36)

To explain the qualitative differences between LHR and B2AR trafficking observed in [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF], we tested two hypotheses here.

• (H1) In the first hypothesis (Figure 7), we modify only the internalization rate α > 0 between LHR and B2AR, with a higher internalization rate for the B2AR, keeping all the remaining parameters the same. This hypothesis is in line with [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF]. LHR and B2AR affect only the endocytosis rate as follows (refer to Eq. ( 35)):

Parameters K 0 α LHR 5 × 10 -1 8 × 10 -5 B2AR 5 × 10 -1 3 × 10 -4 (37)
• (H2) In the second hypothesis (Figure 8), both the internalization rate and the coagulation rate are higher for the B2AR compared to the LHR (H2). LHR and B2AR affect the internalization and the coagulation rate as follows:

Parameters K 0 α LHR 5 × 10 -3 8 × 10 -5 B2AR 5 × 10 -1 3 × 10 -4 (38)
In [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF], mean and variance of the size of the endosomal populations are measured at different time after ligand stimulation, as well as the internalization ratio, which corresponds to the fraction of internalised receptors among the total number of receptors initially present at the plasma membrane. Consistently, we define in the model the internalization ratio as ∥f ∥ 0,1 (t) M0

. In both hypotheses the internalization ratio is well reproduced by the model, compared to [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF]. The major discrepancy of the first hypothesis with the experiments is the variance in size of the endosomal population, which seems too high compared to experimentation. Furthermore the production of large endosomes with LHR seems too high and similar to the ones of B2AR (Figure 7). In the second hypothesis, however, these two discrepancies are not present anymore (Figure 8). These results indicate that the differences in endosome dynamics between LHR and B2AR signaling pathways seem not only due to a difference in internalization rate, but probably also to the coagulation dynamics inside cells. Difference in coagulation dynamics may be explained by differences in endosomes sorting and/or differences in molecular composition of endosomes, which are believed to be of different nature between LHR and B2AR vesicles [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF]. With both hypotheses, we observe a much more efficient cAMP production at the plasma membrane with LA-PTH and a much more efficient cAMP production in the endosomes with PTH 7D (Figures 9 and10).

Consistently with the observation in [START_REF] White | Spatial bias in cAMP generation determines biological responses to PTH type 1 receptor activation[END_REF], both total responses have similar magnitude for the time period of the numerical simulation. However, from the numerical simulation presented in figures 9 and 10, the cAMP production has already reached a "stable" state at T = 20 for LA-PTH, while it keeps increaseasing for PTH 7D. Therefore, a longer time measurement could discriminate between both ligands. Also we could notice a fine kinetic difference between the responses induces by LA-PTH and PTH 7D with the two hypotheses. Indeed, PTH 7D leads to a convex kinetic production of cAMP during the early dynamics, which switches to a concave kinetic at later time. Whereas with LA-PTH the production stays concave all time long. Of course this behaviour may be quite complicated to observe experimentally due to the accuracy of the measures. only by the production rate). In both panels, the dotted blue line shows the endosomal quantity of cAMP, given by ∥f ∥ 0,1 (t), the dashed blue line shows the cAMP quantity at the plasma membrane, given by M (t), and the dashed-dotted purple line shows the total cAMP quantity (∥f ∥ 0,1 (t) + M (t)).

  r, a) = 0.5 • [N r (1.5, 0.15) × N a (0.5, 0.3) + N r (0.5, 0.3) × N a (1.5, 0.15)] , κ (r, r ′ ) = 0.5 .

f 0

 0 (r, a) = 0.5 • [N r (1.5, 0.15) × N a (0.5, 0.3) + N r (0.5, 0.3) × N a (1.5, 0.15)] , κ (r, r ′ ) = 0.5 .

f 0

 0 (r, a) = 0.5 • [N r (1.5, 0.15) × N a (0.5, 0.3) + N r (0.5, 0.3) × N a (1.5, 0.15)] , M 0 = 20 , κ (r, r ′ ) = 0.5 , α (r, a, M ) = 0.1 • M • [N r (0.6, 0.01) × N a (0.3, 0.05) + N r (0.3, 0.05) × N a (0.6, 0.01)] , γ (r, a) = 20 (r -5) 4 1 r>5 + 10 -5 , λ (r, a) = 10 -2 • P r (10, 0) , V (r, a) = 0 , J M (M ) = 0 ,

Figure 4 .

 4 Figure 4. Illustration of the convergence of the numerical scheme (21) in the case of the pure coagulation (case 1, numerical and model parameters given by Eqs. (30)-(31)) . We plot the error as a function of the grid size h, in log-log scale, using the six different norms ∥E(h m , t n )∥ k,l defined in Eq. (28) for 0 ≤ k + l ≤ 2 (see inserted legend).

Figure 5 .

 5 Figure 5. Illustration of the convergence of the numerical scheme (19) in case 2 (numerical and model parameters given by Eqs. (30)-(32)). See legend of Figure 4.

Figure 6 .

 6 Figure 6. Size-distribution of the population density g of endosomal compartement for model (M1) at different times (see inserted legend) and parameters given by Eqs. (33)-[START_REF] Genz | Remarks on algorithm 006: An adaptive algorithm for numerical integration over an Ndimensional rectangular region[END_REF]. On the left hand-side the condition Eq. (6) is satisfied (γ (r) = √ 3.1) whereas on the right hand-side the condition Eq. (6) is not satisfied (γ (r) = 0.7). In both cases, the curves for t = 10 and t = 50 are superimposed.

Figure 7 ..Figure 8 .•

 78 Figure 7. Time evolution of the endosomal size (left panel) and internalization ratio (right panel).We numerically simulate the 2D model given by Eqs. (M2a)-(M2b) with parameters given by. (35)-(36)-(37) (LHR and B2AR affect only the internalization rate). On the left panel, the dashed blue (resp. dotted purple) line represents the mean endosomal size for B2AR (resp. LHR), that corresponds to ∥f ∥ 1,0 (t)/ ∥f ∥ 0,0 (t), and the shaded light area represents its standard deviation. On the right panel, the dashed blue (resp. dotted purple) line represents the internalization ratio for B2AR (resp. LHR), that corresponds to ∥f ∥ 0,1 (t) M 0 .

Figure 9 .

 9 Figure 9. Time evolution of the cAMP production for the LA-PTH ligand (left panel) and the PTH-7D ligand (right panel). We numerically simulate the 2D model given by Eqs. (M2a)-(M2b) with parameters given by. (39)-(40)-(41) (cAMP kinetic parameters differ between PTH 7D and LA-PTH only by the production rate). In both panels, the dotted blue line shows the endosomal quantity of cAMP, given by ∥f ∥ 0,1 (t), the dashed blue line shows the cAMP quantity at the plasma membrane, given by M (t), and the dashed-dotted purple line shows the total cAMP quantity (∥f ∥ 0,1 (t) + M (t)).

Figure 10 .

 10 Figure 10. Time evolution of the cAMP production for the LA-PTH ligand (left panel) and the PTH-7D ligand (right panel). We numerically simulate the 2D model given by Eqs. (M2a)-(M2b) with parameters given by. (39)-(40)-(42) (cAMP kinetic parameters differ between PTH 7D and LA-PTH both by the production rate and capacity of production.). See legend of Figure 9 for details.

https://github.com/JuliaMath/HCubature.jl

-r). We don't consider reactions, e.g. V and J M are taken as null functions. For the numerical scheme, we take

* This work is supported by INRAE Metaprogramme DIGIT-BIO (Digital biology to explore and predict living organisms in their environment) * * This work is supported by a public grant overseen by the French National research Agency (ANR) as part of the ≪ Investissements d'Avenir ≫ program, through the "ADI 2020" project funded by the IDEX Paris-Saclay, ANR-11-IDEX-0003-02 * * * This work is supported by the Italian Ministry of Education, Universities and Research through the MIUR grant "Dipartimenti di Eccellenza 2018-2022" (Project no. E11G18000350001) * * * * This work is partially supported by the spanish State Research Agency and European Regional Development Fund (grant PID2022-137228OB-I00), by Junta de Andalucía and ERDF (grant C-EXP-265-UGR23) and by Modeling Nature Research Unit (project QUAL21-011)

this section is to provide numerical simulations of our model (M2a)-(M2b) that can reproduce these two distinct scenarios.

In order to compare qualitatively our model with the experimental results presented in [START_REF] Jean-Alphonse | Spatially Restricted G Protein-coupled Receptor Activity via Divergent Endocytic Compartments[END_REF], we choose the following parametrization (we recall definitions in Eqs. ( 24)-( 29)):

which we interpret as follows. Endosomes fuse at a constant coagulation kernel κ, whose rate will depend on the receptor. Endosomes are created with a size following a Gaussian law and with a quantity of reactant proportional to their size, and the rate α > 0 will depend on the receptor. Small endosomes are degraded at constant (low) rate, and degradation rates quickly increases for endosomes large enough. Endosomes recycling increase with the surface of endosomes and decreases with their volume ( λ ∝ r

Second effector signaling

The second application of our model we present here concerns the efficacy of second messenger molecules production as a function of the localization of the active receptor. We focus on the production of cAMP induced by the activation of the PTHR. Recent discoveries find that the PTHR may engage cAMP signaling not only at the cell's plasma membrane but also in early endosomes after receptor internalization through endocytosis [START_REF] Vilardaga | Endosomal generation of cAMP in GPCR signaling[END_REF]. Furthermore, the full-length parathyroid hormone (LA-PTH) induces through the activation of PTHR an augmentation of production of cAMP in the endosomes, whereas the ligand PTH 7D (PTHR peptide ligand through amino acid epimerization at position 7 of P T H 1-34 ) induces the production of cAMP at the plasma membrane. In White's paper [START_REF] White | Spatial bias in cAMP generation determines biological responses to PTH type 1 receptor activation[END_REF], the authors show that even if the production place is different, the total amount of cAMP stays the same after some time, a phenomenon that could be named location-biased, and that can have implications for the cellular response. The objective of this section is to provide numerical simulations of our model (M2a)-(M2b) that can reproduce these two distinct scenarios. We could reproduce qualitatively these observations with our model following this parametrization:

which we interpret as follows. Endosomes fuse via a constant coagulation kernel, and are created with a size following a Gaussian law and with a quantity of reactant, independently of the size, taken as a positive fraction of the quantity M present at the plasma membrane. Endosomes recycling and degradation increase with the surface of endosomes and decrease with their volume (λ and γ ∝ r 2 3 -r). cAMP is produced at the plasma membrane at constant rate and linearly degraded, where both the saturation M , and the production rate v M may depend on the ligand. cAMP is produced in endosomes at two different rates, v s , for the endosomes smaller than r and v l , for larger endosomes. Also the amount of saturation depends linearly on the size of the endosomes (caracterized by a portion p of r). Here v s , v l and p may depend on the ligand. For the numerical scheme, we take ∆t = 3 -2 , I r = I a = 30, R = 2000, A = 30 .

(40)

We have two different hypotheses to explain the different qualitative behaviour of the two ligand LA-PTH and PTH 7D described above:

• (H1) Suppose LA-PTH and PTH 7D differs in cAMP production kinetics only in terms of rate (Figure 9), with a higher rate for LA-PTH at the plasma membrane (v P T H7D M < v LA-P T H M ), and a higher rate