Mathematical models for the signaling pathway of the G-protein coupled receptor EP2

Title: Mathematical models for the signaling pathway of the G-protein coupled receptor EP2 Author: Paula Gómez López Advisor: Gemma Huguet Casades Co-advisor: Stefanie Sonner Department: Department of Mathematics Academic year: 2024/2025 Master of Science in Advanced Mathematics and Mathematical Engineering Universitat Polit`ecnica de Catalunya Facultat de Matem`atiques i Estad´ıstica Master in Advanced Mathematics and Mathematical Engineering Master’s thesis Mathematical models for the signaling pathway of the G-protein coupled receptor EP2 Paula G´omez L´opez Supervised by Dr. Gemma Huguet Casades and Dr. Stefanie Sonner January, 2025 This master’s thesis was carried out with the support of the Erasmus+ Internship program from September to December 2024 at Radboud University, under the supervision of Dr. Stefanie Sonner. First of all, I would like to thank Stefanie Sonner for guiding me through this thesis and for answering all my questions. I deeply appreciate Mariya Ptashnyk for all the online meetings and her help with the code. Thanks to my parents for always being there, to Daphne for all the support, and to my friends—both old and new—for making this experience so much better. Abstract This thesis develops mathematical models to study the EP2 signaling pathway, a receptor belonging to the G-protein-coupled receptor (GPCR) family involved in cellular communication. By combining mathematical techniques and biological insights, we aim to better understand the dynamics of EP2 signaling and explore whether spatial effects play an important role. Key contributions include the development of an ODE model to describe the dynamics of cAMP production in EP2 signaling and its extension to incorporate ligand-receptor interactions. Spatial effects are subsequently investigated through reaction-diffusion models. Finally, a novel class of ligand-receptor-based Turing models is introduced to explore receptor clustering via spatial pattern formation, providing a theoretical framework applicable to various receptor systems. This work highlights the power of mathematical modeling in biological signaling and offers tools that extend beyond the EP2 receptor, interconnecting mathematics and biology. Keywords Mathematical modeling, G-protein coupled receptor (GPCR), EP2 signaling, ligand-receptor dynamics, Turing patterns, Schnakenberg type kinetics. 1 Contents 1 Introduction 3 1.1 Biological background ..................................... 4 1.2 Experimental data ....................................... 6 2 Modeling the signaling pathway of EP2 8 2.1 G-protein cycle ......................................... 8 2.2 cAMP production ....................................... 12 2.3 ODE model for cAMP production in EP2 receptor signaling ................. 14 2.3.1 Simulations ...................................... 14 2.4 Modeling ligand-receptor dynamics .............................. 16 2.4.1 Simulations ...................................... 18 3 Spatially heterogeneous model for EP2 signaling 20 3.1 1D reaction-diffusion model in a cross-section of the cell .................. 21 3.1.1 Simulations ...................................... 22 3.2 1D reaction-diffusion model on the cell membrane ...................... 27 3.2.1 Simulations ...................................... 29 4 Ligand-receptor based Turing models 32 4.1 General conditions for diffusion-driven instability ....................... 33 4.1.1 Linear stability in the absence of diffusion ...................... 33 4.1.2 Diffusion driven instability .............................. 34 4.2 Ligand receptor model ..................................... 38 4.3 Ligand-receptor based Turing models ............................. 40 4.3.1 Generalized Schnakenberg model ........................... 42 4.3.2 Model with no feedback receptor production ..................... 51 5 Conclusion and future work 54 6 References 55 A Codes 57 2 Moreover, following [12], we include the reversible Reaction 1in our model, whereas [10] treated this reaction as one-way. The other reactions could also be modeled as reversible, depending on the molecular concentrations and reaction rates. In [10], the Reactions 2,3are treated as enzymatic, that is, the activated receptor catalyzes the dissociation of the G-protein, while GAP catalyzes GTP hydrolysis, modeled by Michaelis-Menten kinetics. In contrast, we assume that Reactions 1and 3are non-enzymatic and only Reaction 2follows MichaelisMenten kinetics. Reactions 1and 3are modeled using the Law of Mass Action, which states that the rate of a reaction is proportional to the product of the concentrations of the reactants [18]. Let [·] denote the concentration of a chemical substance. Moreover, GαGDP sβγ is denoted by αsβγ, GαGDP sby Gαsand GαGTP sby α∗ s. Then, the reaction rates are defined as follows: 1. Association/dissociation of Gαsβγ: The reversible reaction, representing the binding and unbinding of GαGDP sand Gβγ, is not enzyme-catalyzed and is given by: V1=k3[αs][βγ]−k4[αsβγ], where k3and k4are the Gαsβγ association and dissociation rate respectively. 2. Activation of Gαs: The enzymatic reaction, where the activated receptor, [EP2∗], catalyzes the dissociation of Gαsβγ into GαGTP sand Gβγ, is modeled using Michaelis-Menten kinetics. In [25], V2 was initially modeled as: V2=k1[EP2∗][αsβγ] [αsβγ] + K1 , where k1is the Gαsactivation rate and K1is the dissociation constant for EP2 and Gαsβγ. However, based on [10], we simplify the equation to the linear form: V2=k1 K1 [EP2∗][αsβγ]. 3. GTP hydrolysis and deactivation: The deactivation of GαGTP svia GTP hydrolysis, which is not enzymecatalyzed, is given by: V3=k2[α∗ s], where k2is the Gαshydrolysis rate. By Assumption 1, the concentration of activated EP2 receptors, [EP2∗], follows a Hill equation. This means [EP2∗] is determined by the total concentration of EP2 receptors, EP2tot, the dissociation constant, K2for EP2 and PGE2, and the initial concentration of PGE2 added in the experiment. In particular, [EP2∗] varies based on the PGE2 concentration: [EP2∗] = EP2tot[PGE2]n [PGE2]n+K2 . (4) In [12], the Hill coefficient was assumed to be n= 1, corresponding to the classical Michaelis-Menten reaction. This equivalence allows the terms ”Hill equation” and ”Michaelis-Menten reaction” to be used interchangeably when n= 1. In contrast, [25] introduced a higher Hill coefficient, n= 4, to capture the threshold behavior of EP2 for varying ligand concentrations, ensuring the model reflects the qualitative 9 behavior observed in experiments. As shown in Figure 4, the Hill coefficient strongly influences receptor activation, particularly at low PGE2 concentrations. For n= 1, receptor activation saturates quickly as PGE2 increases. For n>1, receptor activation remains minimal at low PGE2 concentrations, creating a threshold effect where weak signals do not trigger a response until a critical concentration is reached. Figure 4: Hill Equation (4) for different Hill coefficients and EP2tot = 0.004µM. Using the mentioned reaction rates, the autonomous ODE system governing the concentrations is: d[αsβγ] dt =V1−V2, (5a) d[βγ] dt =V2−V1, (5b) d[α∗ s] dt =V2−V3, (5c) d[αs] dt =V3−V1, (5d) where as stated before: V1: re-association and dissociation rate of Gαsβγ, V2: dissociation rate of Gαsβγ, V3: hydrolysis rate of GTP on Gαs. The parameter values for the ODE system (5) are taken from [12,25] and are listed in the following table: Parameter Value Units Description k15 s−1Gαsactivation rate k20.07 s−1Gαshydrolysis rate k30.7 µM−1s−1Gαsβγ association rate k418.9 ×10−3s−1Gαsβγ dissociation rate K10.8 µM Dissociation constant for EP2 and Gαsβγ K20.012 µM Dissociation constant for EP2 and PGE2 βγtot 0.005 µM Total Gβγ concentration αtot s2.3 µM Total Gαsconcentration EP2tot 0.004 µM Total EP2 concentration Table 1: Parameters used in the ODE system (5). 10 Regarding mass conservation, two relations can be found: d[αsβγ] dt +d[βγ] dt = 0 =⇒[αsβγ]+[βγ] = βγtot, (6a) d[αsβγ] dt +d[α∗ s] dt =d[αs] dt = 0 =⇒[αsβγ]+[α∗ s]+[αs] = αtot s. (6b) Remark 2.1 (Existence and uniqueness of solutions [21]).Consider the autonomous ODE system: (u′ i(t) = fi(u1, ..., un), i= 1, ..., n, ui(0) = u0 i,i= 1, ..., n. where fi:D→Rnare continuous functions, D⊂Rnis an open set and u(0) = (u1(0), ..., un(0)) ∈D. If the functions fiare continuously differentiable on D, then fiis locally Lipschitz continuous. By the Picard-Lindel¨of theorem, there exists a unique local solution u: [0, T)→Dfor the ODE system, where T>0. Remark 2.2 (Non-negativity criterion).Consider the autonomous ODE system: (u′ i(t) = fi(u1, ..., un), i= 1, ..., n, ui(0) = u0 i≥0, i= 1, ..., n. The solutions ui(t) remain non-negative for all t≥0 if and only if the following condition holds for each i and for all j=i: fi(u1, ..., 0, ..., un)≥0, whenever uj≥0, for all j=i. Remark 2.3 (Existence and uniqueness of solutions of System (5)).The functions on the right-hand sides of the equations are continuously differentiable with respect to [αsβγ], [βγ], [α∗ s], [αs], as they consist of polynomial terms. Therefore, by Remark 2.1, for given initial concentrations, System (5) has a unique local solution around the initial conditions. Remark 2.4 (Non-negativity and boundedness of solutions of System (5)).If [αsβγ] = 0 and [βγ], [α∗ s], [αs]≥ 0, the system satisfies d[αsβγ] dt ≥0. The same holds for the other ODEs. Thus, by the non-negativity criterion in Remark 2.2, all concentrations [αsβγ], [βγ], [α∗ s], [αs] remain non-negative for all t≥0. Furthermore, the mass conservation relations (6) guarantee that the concentrations remain bounded by βγtot and/or αtot s. Using the mass conservation relations (6), we can write the reaction rates V1,V2,V3in terms of [βγ], [α∗ s]. V1=k3[αs][βγ]−k4[αsβγ] = k3[βγ]([βγ]−[α∗ s]−βγtot +αtot s)−k4(βγtot −[βγ]), V2=k1 K1 [EP2∗][αsβγ] = k1 K1 [EP2∗](βγtot −[βγ]), V3=k2[α∗ s]. The system of four ODEs in (5) can be reduced to the system of two ODEs: d[βγ] dt =V2−V1= (k4+k1 K1 [EP2∗])βγtot −k1 K1 [EP2∗] + k4+k3αtot s−k3βγtot [βγ] + k3[βγ][α∗ s]−k3[βγ]2, d[α∗ s] dt =V2−V3=k1 K1 [EP2∗]βγtot −k1 K1 [EP2∗][βγ]−k2[α∗ s]. After solving the system for the two variables [βγ] and [α∗ s], the concentrations of [αsβγ] and [αs] are determined using the mass conservation relations (6). 11 2.2 cAMP production To model cAMP production, we account for both the synthesis by the enzyme AC and the degradation by the PDE enzyme, as mentioned in Section 1. We consider two different equations taken from [25] and [7]. Equation from [25] To capture the positive influence of AC on cAMP levels, we assume that one dominant AC isoform is responsible for cAMP production. A fraction of the total AC, ACtot, is activated by binding to Gα∗ s. We model it using Hill’s equation with n= 1: [AC∗] = ACtot[α∗ s] [α∗ s] + K4 . Some fraction of [AC∗] binds further to Gβγ, forming a ”super-activated” AC∗-βγ complex with superactivation factor C1: [AC∗-βγ]=[AC∗]C1[βγ] [βγ] + K6 . The unbound fraction of [AC∗] is given by: [AC∗ unbound] = [AC∗]−[AC∗-βγ]=[AC∗]K6 [βγ] + K6 . To model the cAMP degradation via PDE enzymes, we use Hill’s equation with n= 1 to describe the concentration of the PDE-cAMP complex: [PDE-cAMP] = PDEtot[cAMP] [cAMP] + Kd , where we focus on the PDE3 and PDE4 isoforms, as they play key roles in regulating cAMP levels [12]. Thus, the equation describing the cAMP concentration, slightly modified from [12] and as given in [25] is: d[cAMP] dt =k9 ACtot[α∗ s] [α∗ s] + K4 C1[βγ] + K6 [βγ] + K6 −k11 PDE4tot[cAMP] [cAMP] + K7 −k12 PDE3tot[cAMP] [cAMP] + K8 . (7) The parameter values are given in Table 2. Equation from [7] Unlike Equation (7), the super-activation term with factor C1is not included here. As a result, the production term for cAMP is given by: Atot[α∗ s] [α∗ s] + Kas . The degradation of cAMP, catalyzed by the PDE enzyme, is modeled using a Michaelis-Menten reaction term: Ptot[cAMP] [cAMP] + Kw . 12 Parameter Value Units Description k96.713 s−1Active cAMP production rate k11 0.72 s−1cAMP degradation rate by PDE4 k12 6.85 s−1cAMP degradation rate by PDE3 K40.2 µM Dissociation constant for AC and Gαs K60.09 µM Dissociation constant for AC and βγ K72.6 µM Dissociation constant for PDE4 and cAMP K80.15 µM Dissociation constant for PDE3 and cAMP ACtot 0.029 µM Total AC concentration PDE4tot 0.115 µM Total PDE4 concentration PDE3tot 0.0025 µM Total PDE3 concentration C111 βγ super-activation factor Table 2: Parameter values for cAMP production in Equation (7), taken from [12]. Hence, the cAMP production is described by: d[cAMP] dt =kw Atot[α∗ s] [α∗ s] + Kas −dw Ptot[cAMP] [cAMP] + Kw . (8) Here, we only consider the main enzymes AC and PDE, which are responsible for cAMP production and degradation, respectively. For simplicity, we assume that a single dominant isoform exists for each enzyme. This assumption is made because we only have qualitative experimental data and no detailed information about the specific enzymes and isoforms, which would be necessary for a more refined modeling approach. The parameter values are given in the following table: Note that since relative concentrations are Parameter Value Units Description kw6.713 s−1Active cAMP production rate Kas 0.2 µM Dissociation constant for AC6 and Gαs dw8.66 s−1cAMP degradation rate by PDE4 Kw1.21 µM Dissociation constant for PDE4 and cAMP Atot 0.0497 µM Total AC concentration Ptot 0.039 µM Total PDE concentration Table 3: Parameter values for cAMP production in Equation (8), taken from [7]. measured in the experiments, we assume that there is no basal cAMP production rate in Equations (7) and (8). 13 2.3 ODE model for cAMP production in EP2 receptor signaling Combining the models for the G-protein activation cycle and cAMP production, we obtain: d[βγ] dt = (k4+k1 K1 [EP2∗])βγtot −k1 K1 [EP2∗] + k4+k3αtot s−k3βγtot[βγ] + k3[βγ][α∗ s]−k3[βγ]2, (9a) d[α∗ s] dt =k1 K1 [EP2∗]βγtot −K1[EP2∗][βγ]−k2[α∗ s], (9b) d[cAMP] dt =kw [α∗ s]Atot [α∗ s] + Kas −dw Ptot[cAMP] [cAMP] + Kw . (9c) where [EP2∗] is given, as in (4), by [EP2∗] = EP2tot[PGE2]n K3+[PGE2]n,n= 4, and the parameters are given in Tables 1and 3. Remark 2.5 (Existence and uniqueness, non-negativity and boundedness of solutions of System (9)).The equations for [βγ] and [α∗ s] are uncoupled from the one for [cAMP]. By Remarks 2.3 and 2.4, Equations (9a) and (9b) have a unique local solution that is non-negative and bounded if the initial conditions satisfy 0≤[βγ]0≤βγtot, 0 ≤[α∗ s]0≤αtot s. Since the reaction terms in the ODE for [cAMP] are continuously differentiable (Kas and Kware positive), Remark 2.1 guarantees the existence of a unique local solution for Equation (9c). Furthermore, by Remark 2.2, the cAMP concentration remains non-negative. However, the solution for [cAMP] is not necessarily bounded. 2.3.1 Simulations Initial conditions The initial conditions for [βγ], [α∗ s] are set to their basal steady-state concentrations. Following [12], we compute the steady-states of the system in the absence of the ligand. To determine these steady-states, we set [PGE2] = 0, which using Equation (4), leads to [EP2∗] = 0. At the steady-state, the reaction rates balance, so V1=V2=V3. Under these conditions we find [α∗ s]0= 0, and for [βγ] we obtain a quadratic equation: k3[βγ]2+ (k4+k3αtot s−k3βγtot )[βγ]−k4βγtot = 0 This equation gives a single biologically relevant positive solution: [βγ]0= 5.81 ×10−5≈6×10−5. Furthermore, we assume an initial concentration of zero for cAMP, [cAMP]0= 0. For the numerical simulations, we use the odeint function in Python, see Appendix A. Figure 5c and 5d show that equations (7) and (8) have different steady-state concentrations for high initial conditions of [PGE2]. However, both equations exhibit similar qualitative behavior, particularly in modeling the cAMP production threshold between low and high ligand concentrations. Given the lack of detailed experimental data on specific enzymes and isoforms, we use Equation (8) for its reduced complexity. The simplifying assumptions in Equation (8) include a single dominant isoform, the omission of the super-activation term and fewer parameters. These assumptions make the model easier to analyze both analytically and numerically, while still capturing the essential threshold behavior for cAMP production. 14 (a) Gβγ concentration. (b) Gαsconcentration. (c) cAMP production using Equation (7). (d) cAMP production using Equation (8). Figure 5: Solutions of System (9). In [25], the threshold between low and high PGE2 concentrations was modeled by increasing the Hill coefficient to n= 4 in Equation (4), which gives the concentration of activated receptors [EP2∗]. A Hill coefficient greater than 1 indicates positive cooperativity, typically suggesting that the receptor has multiple ligand-binding sites. However, GPCRs generally have only a single binding site [25]. This cooperativity could potentially arise from receptor dimerization. While [25] has a more detailed discussion on this threshold mechanism and explores possible mechanisms for this behavior, no definitive biological argument can be given for the Hill coefficient n= 4. Further experiments are needed to determine which Hill exponent best fits the experimental data for critical ligand concentrations near the threshold, specifically [PGE2] ∈0.1, 1. Thus, System (9) successfully reproduces the qualitative experimental dynamics shown in Figure 3. At low ligand concentrations, the EP2 receptor remains inactive, resulting in no cAMP production. However, at higher ligand concentrations, cAMP production increases and stabilizes at a steady-state level. 15 2.4 Modeling ligand-receptor dynamics In [25] and the previous section, the ligand-receptor interaction is assumed to be rapid allowing the use of a quasi-steady-state approximation for the activated receptor concentration (Equation (4)). Here, we aim to extend the model and include the full dynamics of the ligand-receptor interaction to evaluate the validity of this assumption. Our approach is based on the model presented in [9], which mathematically formulates ligand-receptor interaction-driven cell migration in the presence of decoy receptors. By incorporating the ligand-receptor dynamics into System (9), we simulate the full model and compare the original model (9), which uses the steady-state approximation for [EP2∗], with the reformulated model that explicitly includes the ligandreceptor dynamics. According to [7], EP2 receptors do not internalize, so we set the internalization rate ki= 0 for our EP2 receptors, see [9]. The activation of the EP2 receptor, represented by EP2∗, occurs through ligand-receptor binding with the ligand PGE2. This process is described by the reversible reaction: EP2 + PGE2 ka ⇌ kd EP2∗, (10) where kais the association or binding rate constant, and kdis the dissociation or degradation rate constant. Using the law of mass action, previously introduced in Section 2, we derive the following system of ODEs to model the dynamics of each concentration in this reaction. Again, let [·] represent the concentration of each substance. Then, the modified equations from [9], governing the receptor-ligand dynamics are given by: d[PGE2] dt =U1−U2=kd[EP2∗]−ka[PGE2][EP2], (11a) d[EP2∗] dt =U2−U1=ka[PGE2][EP2] −kd[EP2∗], (11b) d[EP2] dt =U1−U2=kd[EP2∗]−ka[PGE2][EP2]. (11c) First, we outline the mass conservation relations of the system d[EP2∗] dt +d[EP2] dt = 0 =⇒[EP2∗] + [EP2] = EP2tot, (12a) d[PGE2] dt +d[EP2∗] dt = 0 =⇒[PGE2] + [EP2∗] = PGE2tot. (12b) Thus, the system of three dependent equations can be reduced to one independent equation, describing the dynamics of the whole System (11). However, we only make use of the activated and non-activated receptor mass conservation (12a). 16 Secondly, we look for the steady-states of the system. The receptor-ligand system reaches the steadystate when the time derivatives in Equations (11) become zero, that is, the rates U1and U2become equal. U2−U1= 0 ⇐⇒ ka[PGE2][EP2]−kd[EP2∗] = 0 ⇐⇒ ka[PGE2](EP2tot −[EP2∗])−kd[EP2∗]=0 ⇐⇒ ⇐⇒ (ka[PGE2] + kd)[EP2∗] = kaEP2tot[PGE2] ⇐⇒ [EP2∗] = kaEP2tot[PGE2] ka[PGE2] + kd =EP2tot[PGE2] [PGE2] + kd/ka . Hence, the steady-state expression for the activated EP2 receptor is given by: [EP2∗] = EP2tot[PGE2] [PGE2] + K, where K=kd ka. This expression for [EP2∗] matches Equation (4), with n= 1. To generalize the system for any Hill coefficient n, particularly n= 4, to model the threshold behavior, we replace [PGE2] by [PGE2]n in Equations (11). This results in the desired system for receptor-ligand dynamics: d[PGE2] dt =U1−U2=kd[EP2∗]−ka[PGE2]n[EP2], (13a) d[EP2∗] dt =U2−U1=ka[PGE2]n[EP2] −kd[EP2∗], (13b) d[EP2] dt =U1−U2=kd[EP2∗]−ka[PGE2]n[EP2]. (13c) To estimate the association constant ka, we use the dissociation rate kd= 0.0058s−1from [9] and the equilibrium constant K2=K= 0.012µMfrom [25], and we compute ka=kd/K2. The final system of equations for the simulation, taking into account ligand-receptor dynamics, includes the dynamics of [βγ], [α∗ s] and [cAMP]: d[PGE2] dt =kd[EP2∗]−ka[PGE2]n[EP2], (14a) d[EP2∗] dt =ka[PGE2]n[EP2] −kd[EP2∗], (14b) d[EP2] dt =kd[EP2∗]−ka[PGE2]n[EP2], (14c) d[βγ] dt = ( k1 K1 [EP2∗] + k4)βγtot −k1 K1 [EP2∗] + k4+k3αtot s−k3βγtot[βγ] + k3[βγ][α∗ s]−k3[βγ]2, (14d) d[α∗ s] dt =k1 K1 [EP2∗]βγtot −k1 K1 [EP2∗][βγ]−k2[α∗ s], (14e) d[cAMP] dt =kw [α∗ s]Atot [α∗ s] + Kas −dw Ptot[cAMP] [cAMP] + Kw . (14f) Remark 2.6 (Existence and uniqueness, non-negativity and boundedness of System (14)).By Remarks 2.1 and 2.2, Equations (14a), (14b) and (14c) have unique local solutions that are non-negative and bounded by EP2tot and/or PGE2tot for non-negative initial conditions. Similarly, by Remark (2.5), Equations (14d), (14e) and (14f) also admit unique, non-negative local solutions. 17 2.4.1 Simulations Initial conditions The initial conditions for [βγ], [α∗ s] and [cAMP] are the same as in System (9): [α∗ s]0= [cAMP]0= 0 and [βγ]0= 6 ×10−5. We consider four different initial values for the ligand concentration: [PGE2]0∈ {0.01, 0.1, 1, 10}. The initial concentration of activated receptors is set to zero, while the non-activated receptor concentration is determined using the mass conservation relation (12). Thus, [EP2∗]0= 0 and [EP2]0= EP2tot. The plots in Figure 6confirm the validity of the rapid ligand-receptor dynamics assumption, justifying the use of the quasi-steady-state approximation for [EP2∗]. The EP2 receptors activate quickly, resulting in no qualitative differences between the plots in Figures 5and 6. Consequently, System (9) is sufficient to describe the EP2 signaling pathway. Moreover, the ligand concentration [PGE2] remains nearly constant throughout the process. This stability arises because the binding of PGE2 to EP2 involves only a small fraction of the total ligand, leading to negligible changes in [PGE2] during receptor interaction. While the equations for [PGE2] and [EP2] are structurally identical (see Equations (14a) and (14c)), their initial conditions differ significantly. Specifically, [EP2] starts at 0.004, while [PGE2] takes values of 0.01, 0.1, 1, and 10. This difference in initial concentrations significantly influences the dynamics observed in Figure 6. This highlights how initial conditions can determine the behavior. We remark here that the threshold observed in [25] occurs due to the nonlinearity introduced by the ligand concentration in the steady-state approximation of the activated receptors, EP2∗. Mathematically, this is characterized by the dependence of [EP2∗] on [PGE2]n. For small ligand concentrations, the term [PGE2]nbecomes negligible, preventing receptor activation. Conversely, for sufficiently large ligand concentrations ([PGE2] ≫1), the term [PGE2]ndominates, rapidly saturating receptor activation and driving [EP2∗] to its steady-state value. 18 b(k)=           u(k),n 0+ ∆t·f(k),n 0+∆t ∆xg(k),n+1 0 u(k),n j+ ∆t·f(k),n j . . . u(k),n Nx−1+ ∆t·f(k),n Nx−1+∆t ∆xg(k),n+1 L           . We then solve the linear system A(k)u(k),n+1 =b(k)at each time step using spsolve in Python. Note that the terms g(k),n+1 0and g(k),n+1 Lcorrespond to the fluxes at the boundaries, which in System (15), depend on the solution of the uncoupled ligand-receptor ODE system. These boundary fluxes must be updated at each time step using the solution of the ODE system. Then, this ODE system must be solved at each time step using the Backward Euler implicit method: u(k),n+1 =u(k),n+ ∆t·f(k)(tn+1,u(k),n+1). The solutions of this ODE system at each time step are used to update the flux boundary conditions for the reaction-diffusion system. The original ODE model, designed to match the experimental data, accurately represents the dynamics of EP2 signaling. Adding intracellular diffusion to the model, while keeping the same parameter values and equations as in the original ODE model, shows that the key signaling dynamics remain localized at the membrane. Figures 8,9and 10 as well as System (15) show cAMP production occurs on the cell membrane, with its concentration decaying in the intracellular space. This decay, however, is not significantly influenced by diffusion. This suggests that intracellular diffusion is negligible, and the simpler ODE model is sufficient to describe the system. Additionally, the membrane dynamics closely resemble those described by the original ODE model. That is, System 15 preserves the switch-like behavior for different initial conditions. Particularly, the small production of cAMP observed in Figure 8reflects boundary conditions influencing [βγ] and [α∗ s], preserving the switch-like behavior under different initial ligand conditions. However, when adding initial ligand concentrations of 1µM or 10µM, the production of cAMP on the cell membrane increases, see Figures 9and 10. Given that most concentrations are localized at the boundaries of the domain, this motivated the development of a second PDE model focused on the 1-dimensional cell membrane to better explore spatial effects in that region. 25 Figure 8: Left: initial conditions for system (15). Right: Solution of (15) for [PGE2]0= 0.1µM at time T=100. Figure 9: Solution of (15) for [PGE2]0= 1µM at time T=100. Figure 10: Solution of (15) for [PGE2]0= 10µM at time T=100. 26 3.2 1D reaction-diffusion model on the cell membrane The reaction-diffusion model in the previous section showed that diffusion of substrates within the intracellular space is negligible. This suggests that the relevant signaling processes are localized at the cell membrane, where the relevant substrates are concentrated. Consequently, we introduce another reaction-diffusion model on the simplified 1-dimensional domain representing the cell membrane, denoted by Ωc= (0, L). This approach is justified by the membrane localization of EP2 receptors and adenylyl cyclase (AC), the enzyme responsible for cAMP synthesis [8,13]. To incorporate spatial diffusion along the membrane, diffusion terms are introduced for each substrate in the system, including ligands, receptors (both activated and non-activated forms), and G-protein subunits. Smaller molecules, such as cAMP and ligands, are assumed to diffuse faster than larger, membrane-bound structures like receptors. The diffusion coefficients are the same as the ones used previously. Hence, diffusion occurs along the cell membrane, allowing substrates to move within the domain Ωcand interact through the defined reaction kinetics. All reactions are assumed to happen in the domain. That is, Reactions (1), (2) and (3) corresponding to the G-protein cycle, ligand-receptor binding (10) and cAMP production and degradation by enzymes AC and PDE respectively, are considered in the domain Ωc. Regarding the boundary conditions, we apply homogeneous Neumann boundary conditions on the boundary Γc={0, L}. This choice represents the assumption of no flux of chemicals across the boundary, meaning that the molecules (ligand, receptors and G-protein subunits) do not enter or leave the defined 1-dimensional membrane domain, Ωc. These boundary conditions are reasonable because they reflect the restriction of signaling events and molecules to the cell membrane, where the main processes of EP2 signaling occur. Alternatively, periodic boundary conditions can be applied, to emphasize the closed-curve structure of the cell membrane. Figure 11: 1-dimensional cell-membrane. 27 We take the ODE System (14) and add diffusion for every substrate. This model is represented by the following system of PDEs over the domain, Ωc: ∂[PGE2] ∂t=D1 ∂2[PGE2] ∂x2+U1−U2=D1 ∂2[PGE2] ∂x2+kd[EP2∗]−ka[PGE2]n[EP2], (20a) ∂[EP2∗] ∂t=D3 ∂2[EP2∗] ∂x2+U2−U1=D3 ∂2[EP2∗] ∂x2+ka[PGE2]n[EP2] −kd[EP2∗], (20b) ∂[EP2] ∂t=D3 ∂2[EP2] ∂x2+U1−U2=D3 ∂2[EP2] ∂x2+kd[EP2∗]−ka[PGE2]n[EP2], (20c) ∂[βγ] ∂t=D2 ∂2[βγ] ∂x2+V2−V1=D2 ∂2[βγ] ∂x2+k4βγtot +k1 K1 βγtot[EP2∗] −(k4+k3αtot s−k3βγtot)[βγ]−k1 K1 [EP2∗][βγ] + k3[βγ][α∗ s]−k3[βγ]2, (20d) ∂[α∗ s] ∂t=D2 ∂2[α∗ s] ∂x2+V2−V3=D2 ∂2[α∗ s] ∂x2+k1 K1 βγtot [EP2∗]−k1 K1 [EP2∗][βγ]−k2[α∗ s], (20e) ∂[cAMP] ∂t=D1 ∂2[cAMP] ∂x2+kw [α∗ s]Atot [α∗ s] + Kas −dw Ptot[cAMP] [cAMP] + Kw , (20f) with homogeneous Neumann boundary conditions on the boundary Γc={0, L}. As in the previous section, we introduce kfor compactness of notation: k=(1, on Γ1 c={x= 0}, 0, on Γ2 c={x=L}. Then, the boundary conditions on Γcare given by: (−1)kD1 ∂[PGE2] ∂x= 0, (−1)kD3 ∂[EP2∗] ∂x= 0, (−1)kD3 ∂[EP2] ∂x= 0, (−1)kD2 ∂[βγ] ∂x= 0, (−1)kD2 ∂[α∗ s] ∂x= 0, (−1)kD1 ∂[cAMP] ∂x= 0. Remark 3.3.Note that we consider the same diffusion coefficient for the G-proteins, allowing us to apply the mass conservation relations (6) and reduce the system of four equations to two equations. We also assume equal diffusion coefficients for the activated and non-activated receptors. However, we do not apply the mass conservation relation (12) in this case. Remark 3.4 (Existence, uniqueness, non-negativity and boundedness of solutions of System (20)).For the existence and uniqueness of solutions, we refer to [26]. The non-negativity criterion in Remark 2.2 also applies to reaction-diffusion systems (see [26]). Since the reaction terms remain unchanged, the solutions remain non-negative. 28 3.2.1 Simulations Initial Conditions As in the ODE model, we assume basal initial concentrations for all substrates: [βγ](x, 0) = 6×10−5, [α∗ s](x, 0) = 0, [EP2∗](x, 0) = 0, [PGE2](x, 0) ∈ {0.01, 0.1, 1, 10},x∈(0, L). To reflect biological evidence that EP2 receptors are clustered in specific regions along the cell membrane [6], we introduce a spatially heterogeneous initial condition for [EP2] with concentration peaks. We model this distribution as a sum of Gaussian functions centered at distinct points along the domain: [EP2](x, 0) = N X i=1 hiexp −(x−pi)2 2σ2, where Nis the number of peaks, hiis the amplitude of the i-th peak, piis its center position, and σ controls the width of the peaks, representing receptor clustering. To ensure that the results of this PDE model are comparable to those of the ODE model, we take the height of the peaks to be equal to EP2tot = 0.004. This choice of the initial conditions provides a heterogeneous and biologically realistic initial distribution for [EP2]. Moreover, the initial conditions are comparable with the values from the ODE model. Implementation The discretization and implementation of System (20) follows the same approach as in the previous section. In the interior points, the solution is approximated using Equation (16). However, in this case, we impose homogeneous (zero-flux) Neumann boundary conditions, which give: Dk ∂u(k),n+1 j ∂xj=0 ≈Dk u(k),n+1 0−u(k),n+1 −1 ∆x= 0 ⇒u(k),n+1 −1=u(k),n+1 0, Dk ∂u(k),n+1 j ∂xj=Nx−1≈Dk u(k),n+1 Nx−u(k),n+1 Nx−1 ∆x= 0 ⇒u(k),n+1 Nx=u(k),n+1 Nx−1. Using these boundary conditions, we can rewrite the equations for j= 0 and j=Nx−1 as follows: At j= 0 : (1 + λk)u(k),n+1 0−λku(k),n+1 1=u(k),n 0+ ∆t·f(k),n 0. At j=Nx−1 : −λku(k),n+1 Nx−2+ (1 + λk)u(k),n+1 Nx−1=u(k),n Nx−1+ ∆t·f(k),n Nx−1. System (20) is discretized into a linear system of equations: A(k)u(k),n+1 =b(k), where A(k)is the finite difference matrix (19), and b(k)is the modified right-hand side vector with zero-flux conditions, given by: b(k) j=u(k),n j+ ∆t·f(k),n j. 29 When simulating the system with homogeneous initial conditions (including [EP2] = EP2tot = 0.004µM), the model exhibits a clear threshold behavior, like in the ODE system. As shown in Figure 12, for a PGE2 concentration of 1µM or lower, receptor activation does not occur, and there is no cAMP production. However, when [PGE2] reaches 1µM, receptors become activated, and cAMP production reaches a homogeneous steady state along the domain. Figure 12: Left: PGE2 = 0.1µM. Right: PGE2 = 1µM. To explore spatial effects, we consider non-homogeneous initial conditions for the non-activated receptors, inspired by evidence of receptor clustering on the cell membrane [6]. The initial condition models two receptor clusters. Again, Figures 13,14 and 15 show a threshold between low and high concentrations of PGE2. However, numerical simulations show that the diffusion terms homogenize the substrate distributions over time, eliminating receptor clusters, see Figure 14 at T=100 and 15 at T=200. That is, all concentrations tend to have a homogeneous distribution on the domain. As shown in Figures 13,14, and 15, there is a threshold between low and high concentrations of PGE2. However, numerical simulations show that the diffusion terms gradually homogenize the substrate distributions over time, dissipating the receptor clusters. This behavior is shown in Figure 14 at T= 100 and Figure 15 at T= 200, where all concentrations trend toward a uniform distribution across the domain. Figure 13: Left: initial conditions for System (20). Right: concentrations of System (20) for [PGE2]0= 0.1µM. 30 Figure 14: Left: initial conditions. Right: concentrations of System (20) for [PGE2]0= 1µM at T=100. Figure 15: Left: initial conditions. Right: concentrations of System (20) for [PGE2]0= 1µM at T=200. In summary, System (20) reproduces the switch-like behavior observed in the ODE system, capturing the threshold-dependent activation of receptors and cAMP production. However, the system is unable to model EP2 receptor clustering along the cell membrane. In the next section, we introduce ligand-receptorbased Turing models to explore their potential in modeling the receptor clustering on the cell membrane. Therefore, while both reaction-diffusion models (15) and (20) are mathematically valuable, this extension does not provide significant improvements over the ODE model for EP2 signaling. This suggests on the one hand that intracellular diffusion is negligible, making the ODE model (9) sufficient, and on the other hand that the second model, (20), needs to be modified to account for spatial clustering. 31 4. Ligand-receptor based Turing models In the previous section, two reaction-diffusion models were considered to study spatial diffusion of substrates, in the cytoplasm and on the cell membrane. While these models successfully reproduce the threshold behavior observed in experimental data, model (20), proposed on the cell membrane, failed to capture receptor clustering, as experimentally observed by fluorescence microscopy. This difference highlights the limitations of the previous models in explaining spatial patterns on the cell membrane. In this section, we explore a mechanism for producing spatial patterns, focusing on receptor clustering through ligand-receptor interactions. First, we introduce the theory of Turing patterns, which provides the understanding of spatial pattern formation. Next, we justify why the previous models could not generate patterns. Finally, we propose a new approach, proposing a class of models capable of producing patterns under specific simplifying assumptions and conditions. These models have the potential to be applied to EP2 signaling when additional experimental data becomes available and may also find broader applications in other signaling systems. This work is inspired by previous studies of ligand-receptor interactions in different biological systems [11,16,15,17,27,3]. 32 4.1 General conditions for diffusion-driven instability Consider the system of differential equations with zero flux boundary conditions and initial conditions:                ut(x,t)=∆u(x,t) + γfu(x,t), v(x,t),x∈Ω, t>0, vt(x,t) = d∆v(x,t) + γgu(x,t), v(x,t),x∈Ω, t>0, (n· ∇) u v!= 0, on ∂Ω, t>0, u(x, 0) = u0(x), v(x, 0) = v0(x), x∈Ω. (21) where ∂Ω is a closed boundary of the bounded and connected domain Ω ⊂Rnand nis the unit outward normal to ∂Ω. The choice of zero flux boundary conditions is important since these conditions imply no external input. If this was not the case, the boundary conditions on uand vwould affect directly the spatial pattern [18]. We use homogeneous Neumann (zero-flux) boundary conditions because they are mathematically simple to handle, biologically relevant (modeling impermeable boundaries with no flux across the domain), and allow us to focus on the self-organization of patterns without external input. Unlike fixed (Dirichlet) boundary conditions, which can impose spatial patterns directly, zero-flux conditions ensure that any observed spatial structures arise solely from the internal dynamics of the reaction-diffusion system, making them ideal for studying Turing instabilities and emergent patterns. However, Dirichlet boundary conditions would not greatly alter the following analysis [18,11]. Turing instability appears when a reaction-diffusion system has a stable homogeneous steady state in the absence of diffusion, which loses its stability in the presence of diffusion such that spatial patterns emerge. 4.1.1 Linear stability in the absence of diffusion A spatially homogeneous steady-state (u0,v0) of the System (21) satisfies: ut=γf(u,v) = 0, vt=γg(u,v) = 0. We linearize the system around (u0,v0) by introducing the translated function z= (z1,z2)Twith z1=u−u0 and z2=v−v0. Then, the linearized system around the steady-state becomes: zt=γAz, where A=fufv gugv(u0,v0) =fu(u0,v0)fv(u0,v0) gu(u0,v0)gv(u0,v0), is the Jacobian evaluated at the point (u0,v0). From now on, we write the partial derivatives evaluated at the steady-state without their arguments for simplicity. The steady-state of the linearized system is stable, if Re(λ)<0 for all eigenvalues of A, which for a 2-dimensional system is ensured by the conditions: tr(A) = fu+gv<0, det(A) = fugv−fvgu>0. 33 Remark 4.1.The trace and determinant determine the eigenvalues and vice versa: tr(A) = λ1+λ2, det(A) = λ1λ2. 4.1.2 Diffusion driven instability Now, we aim to derive conditions for the system to become unstable under spatial perturbations introduced by diffusion. Consider the full reaction-diffusion system (21). By applying the same perturbation procedure, we linearize the reaction-diffusion system around the equilibrium point z= (0, 0) to obtain the following linearized system: zt=γAz + D∆z=γfu(u0,v0)fv(u0,v0) gu(u0,v0)gv(u0,v0)z+1 0 0d∆z. (22) Using the method of separation of variables, we look for solutions of the form z(x,t) = T(t)X(x), where T(t) is a temporal function and X(x) is a spatial function. Inserting this ansatz into the linearized system (22) we get T′(t)X(x) = γAT(t)X(x)+DT(t)∆X(x). Dividing by T(t)X(x), we obtain: T′(t) T(t)=γAX(x) + D∆X(x) X(x). Here, the left-hand side is a function of tonly, while the right-hand side depends on xonly. The only way these expressions can be equal is if both are constant. Thus, we have: T′(t) T(t)=γAX(x) + D∆X(x) X(x)=˜ λ, for some constant ˜ λ∈R. Hence, the separation of variables leads to two distinct problems. First, the temporal problem is given by the equation T′(t) = ˜ λT(t), (23) which describes the time evolution of the system, and has the general solution: T(t) = T(t0)e˜ λt, where the exponent ˜ λdetermines the temporal growth or decay of the solution. Second, the spatial problem is given by the elliptic eigenvalue problem: (γAX(x) + D∆X(x) = ˜ λX(x), x∈Ω, (n· ∇)X(x) = 0, x∈∂Ω. 34 •Membrane-bound ligand-receptor complex formation is slow compared to the binding and turnover dynamics. Hence, diffusion is not considered in the equation for the ligand-receptor complex C. •Fast dynamics of the receptor-ligand complex, compared to the other chemicals. Hence, we introduce a quasi-steady-state approximation for the complex concentration, −δC[C] + kon[R]m[L]n−koff[C]=0 ⇐⇒ [C] = kon δC+koff [R]m[L]n= Γ[R]m[L]n, where Γ = kon δC+koff . Hence, the concentration of bound/activated receptors, [C], is proportional to [R]m[L]n, Assumption 2: Complex linear dependent receptor upregulation. The rate of ligand-receptor-dependent receptor upregulation depends linearly on the ligand-receptor complex concentration, [C]: µ([C]) = v[C] = vΓ[R]m[L]n. Note that if we expect a saturation of the response for higher ligand-receptor concentrations, we can consider a Hill function, however, this leads to a smaller Turing space [11], µ([C]) = H(µ([C]), K) = H(Γ[R]m[L]n,K). Under these assumptions, the model for the ligand and receptor dynamics reduces to: ˙ [L] = DL∆[L] + ρS−δL[L]−n(kon[R]m[L]n−koff[C]) = =DL∆[L] + ρL−δL[L]−nδC[C] = =DL∆[L] + ρL−δL[L]−nδCΓ[R]m[L]n, ˙ [R] = DR∆[R] + ρR+µ([C]) −δR[R]−m(kon[R]m[L]n−koff[C]) = =DR∆[R] + ρR+vΓ[R]m[L]n−δR[R]−mδC[C] = =DR∆[R] + ρR+−δR[R]+(v−mδC)Γ[R]m[L]n. These equations model the dynamics of ligand and receptor concentrations, incorporating feedback production and receptor-ligand interactions. In [11], it is shown that under specific assumptions, the previous system converges to the Schnakenberg model. These assumptions are as follows: •Receptor-independent ligand degradation is negligible compared to receptor-dependent degradation: δL[L]<< nδcΓ[R]m[Ln=⇒δL= 0. •The stoichiometry of the ligand-receptor interaction in [11] is m= 2, n= 1, corresponding to one receptor dimer binding two monomeric ligands. However, here, we consider general integers m,n, to determine the possible combinations. •In [11], the feedback production linearly depends on the complex concentration with coefficient v= (m+n)δC. However, here, we consider two cases to investigate the impact of feedback production: v=(0, if µ([C]) = 0 (m+n)δC, if µ([C]) = v[R]m[L]n This allows for linear feedback production or no feedback at all. 41 We redefine the variables and parameters to simplify notation. Let U=U(X,τ) be the receptor concentration and V=V(X,τ) the ligand concentration. Moreover, we consider the following rescaled parameters: ρR=k1,δR=k2,k5=vΓ, k3=δCΓ, k4=ρL,δL=k6= 0. With this new notation, the system becomes: ∂U ∂τ =DU∆U+k1−k2U+ (k5−mk3)UmVn, ∂V ∂τ =DV∆V+k4−nk3UmVn. As mentioned, we distinguish the following two cases: •Generalized Schnakenberg model We consider linearly complex dependent receptor upregulation as in [11]: µ([C]) = v[R]m[L]n. Then, v= (m+n)δc, which implies, k5= (n+m)k3. Hence, the system results in a generalized Schnakenberg model: ∂U ∂τ =DU∆U+k1−k2U+nk3UmVn, ∂V ∂τ =DV∆V+k4−nk3UmVn. (28) •No feedback production We do not consider complex dependent receptor production: µ([C]) = 0. Here, v= 0, so, k5= 0. The system becomes: ∂U ∂τ =DU∆U+k1−k2U−mk3UmVn, ∂V ∂τ =DV∆V+k4−nk3UmVn. (29) 4.3.1 Generalized Schnakenberg model First, to have a model independent of the unit system and have a reduced number of parameters, we rewrite System (28) using the following dimensionless variables and parameters: u=Unk3 k21 m+n−1,v=Vnk3 k21 m+n−1,t=DUτ L2,x=X L, d=DV DU ,a=k1 k2nk3 k21 m+n−1,b=k4 k2nk3 k21 m+n−1,γ=L2k2 DU . (30) The left-hand sides of the equations in System (28) remain: ∂U ∂τ =∂ ∂τ "k2 nk31 m+n−1 u#=k2 nk31 m+n−1∂u ∂t ∂t ∂τ =DU L2k2 nk31 m+n−1∂u ∂t, ∂V ∂τ =∂ ∂τ "k2 nk31 m+n−1 v#=k2 nk31 m+n−1∂v ∂t ∂t ∂τ =DU L2k2 nk31 m+n−1∂v ∂t. 42 The Laplacians are given by: DU∆U=DU ∂2U ∂X2=DU ∂ ∂X∂U ∂X=DU ∂ ∂X"k2 nk31 m+n−1∂u ∂x ∂x ∂X#=DU L2k2 nk31 m+n−1∂2u ∂x2, DV∆V=DV ∂2V ∂X2=DV ∂ ∂X∂V ∂X=DV ∂ ∂X"k2 nk31 m+n−1∂v ∂x ∂x ∂X#=DV L2k2 nk31 m+n−1∂2v ∂x2. The reaction terms are as follows: f(U,V) = f k2 nk31 m+n−1 u,k2 nk31 m+n−1 v!=k1−k2k2 nk31 m+n−1 u+nk3k2 nk3m m+n−1k2 nk3n m+n−1 umvn, g(U,V) = g k2 nk31 m+n−1 u,k2 k51 m+n−1 v!=k4−nk3k2 nk3m m+n−1k2 nk3n m+n−1 umvn. Then, we rewrite the equations in the dimensionless form: DU L2k2 nk31 m+n−1∂u ∂t=DU L2k2 nk31 m+n−1∂2u ∂x2+k1−k2k2 nk31 m+n−1 u−nk3k2 nk3m m+n−1k2 nk3n m+n−1 umvn, DU L2k2 nk31 m+n−1∂v ∂t=DV L2k2 nk31 m+n−1∂2v ∂x2+k4−nk3k2 nk3m m+n−1k2 nk3n m+n−1 umvn. Simplifying the terms we obtain: ∂u ∂t=∂2u ∂x2+L2 DUnk3 k21 m+n−1 k1−L2 DU k2u−L2 DU k2umvn, ∂v ∂t=DV DU ∂2v ∂x2+L2 DUnk3 k21 m+n−1 k4−L2 DU k2umvn. Hence, we get the dimensionless system for the ligand-receptor dynamics: ∂u ∂t= ∆u+γ(a−u+umvn), ∂v ∂t=d∆v+γ(b−umvn), (31) where a,b,d,γ > 0. This system generalizes the Schnakenberg model for arbitrary stoichiometric exponents mand n. Note that for m= 2, n= 1, we recover the well-known Schnakenberg model, with reaction terms: f(u,v) = a−u+u2v,g(u,v) = b−u2v. Remark 4.3 (Existence, uniqueness and non-negativity of solutions of System (31)).For the existence and uniqueness of solutions, we refer to [26]. From the generalization of Remark 2.2 (see [26]), solutions of the generalized Schnakenberg model remain non-negative. Pattern formation The goal is to determine the combinations of stoichiometric exponents in the ligand-receptor binding reaction (27) that lead to pattern formation. To do this, we follow Turing’s theory as presented in Section 4.1, performing a linear stability analysis first in the absence of diffusion and then incorporating diffusion-driven 43 instabilities. Homogeneous steady-states Consider System (31) in the absence of diffusion: ut=γf(u,v) = γ(a−u+umvn), vt=γg(u,v) = γ(b−umvn). The homogeneous steady-states (u∗,v∗) are determined by solving: f(u,v) = 0, g(u,v) = 0. The resulting homogeneous steady-states are: (u∗,v∗) = a+b,b (a+b)m1/n!. (32) These homogeneous steady-states, are positive and real for any a,b>0, ensuring their biological relevance in the context of the model. Stability of the steady-states and Turing analysis To analyze the stability of the steady-states, we compute the Jacobian matrix evaluated at (u∗,v∗): J(u∗,v∗) = fu(u,v)fv(u,v) gu(u,v)gv(u,v)(u∗,v∗) =−1 + mum−1vnnumvn−1 −mum−1vn−numvn−1(u∗,v∗) , where: fu(u∗,v∗) = −1 + mb a+b=(m−1)b−a a+b,fv(u∗,v∗) = n(a+b)m nbn−1 n>0, gu(u∗,v∗) = −mb a+b<0, gv(u∗,v∗) = −n(a+b)m nbn−1 n<0. To show stability in the absence of diffusion and diffusion-driven instability, we use the necessary conditions for Turing instabilities derived in Section 4.1.2. 1. fu+gv<0 Substituting the expressions for fuand gvat (u∗,v∗), we have: fu+gv=−1 + mb a+b−n(a+b)m nbn−1 n<0. Rearranging, this gives: mb a+b<1 + n(a+b)m nbn−1 nor equivalently (m−1)b−a<n(a+b)m n+1bn−1 n. 2. fugv−fvgu>0 Expanding the determinant at (u∗,v∗), we find: fugv−fvgu=−−1 + mb a+bn(a+b)m nbn−1 n+mb a+bn(a+b)m nbn−1 n. 44 Simplifying, this reduces to: fugv−fvgu=n(a+b)m nbn−1 n>0, which is always satisfied, since a,b>0. 3. dfu+gv>0 Substituting the derivatives at (u∗,v∗), we obtain: dfu+gv=d−1 + mb a+b−n(a+b)m nbn−1 n>0. Rearranging, this gives: dmb a+b>d+n(a+b)m nbn−1 nor equivalently d(m−1)b−da >n(a+b)m n+1bn−1 n. 4. d= 1 5. fugv<0 Expanding fugvat (u∗,v∗), we find: fugv=−−1 + mb a+bn(a+b)m nbn−1 n<0. This implies: 1−mb a+b<0⇐⇒ 1<mb a+b⇐⇒ a<(m−1)b, requiring m>1, in order to have a>0. 6. (dfu+gv)2−4d(fugv−fvgu)>0 Evaluating at (u∗,v∗), d−1 + mb a+b−nb n−1 n(a+b)m n2 −4dnb n−1 n(a+b)m n>0. Summarizing, for Turing instability to occur, the following conditions must be satisfied: 1. fu+gv<0 =⇒(m−1)b−a<n(a+b)m n+1bn−1 n 2. fugv−fvgu>0 (always satisfied) 3. dfu+gv>0 =⇒d(m−1)b−da >n(a+b)m n+1bn−1 n 4. d= 1 5. fugv<0 =⇒a<(m−1)b=⇒m>1 6. (dfu+gv)2−4d(fugv−fvgu)>0 =⇒d(−1 + mb a+b)−nb n−1 n(a+b)m n2>4dnb n−1 n(a+b)m n 45 To determine which pairs of stoichiometric numbers (m,n) give rise to Turing patterns in System (26), we start by noting that the condition 5implies that no Turing patterns can occur if m= 1. Therefore, we focus on cases where m>1. For example, the pair (m= 2, n= 1), corresponding to the Schnakenberg model, is known to satisfy the necessary conditions for Turing patterns for certain parameter values and to produce such patterns. To systematically determine which other pairs (m,n) allow Turing patterns, we first analyze the necessary conditions in 4.3.1 for ranges of parameters and then perform numerical simulations. Here, we consider all 12 combinations of m∈ {2, 3, 4}and n∈ {1, 2, 3, 4}. For each combination, we check the necessary Turing conditions 4.3.1, over a range of values for the parameters a∈(0, 5), b∈(0, 5) and d∈(5, 100). The results show that all these stoichiometric combinations satisfy the necessary conditions for Turing patterns. Furthermore, numerical simulations confirm the emergence of Turing patterns for each pair, see Figures 19,20 and 21. Now we aim to determine the Turing space for the parameters aand b, while allowing dto vary [4]. To do this, we modify the expressions in 4.3.1 and we write them in terms of aand u∗, given in (32). We just reformulate conditions 1,3and 6from the Turing necessary conditions 4.3.1. The second condition is always satisfied, while conditions 4and 5are derived directly from condition 3. 1. fu+gv<0 =⇒(m−1)(u∗−a)−a<n(u∗)m n(u∗−a)n−1 n⇐⇒ a>u∗ m(m−1) −n(u∗)m n(u∗−a)n−1 n, 3. dfu+gv>0 =⇒d(m−1)u∗−dma >n(u∗)m n+1(u∗−a)n−1 n⇐⇒ a<u∗ m(m−1) −n d(u∗)m n(u∗−a)n−1 n, 6. (dfu+gv)2−4d(fugv−fvgu)>0 =⇒d(m−1) −m u∗a−n(u∗)m n(u∗−a)n−1 n2>4nd(u∗)m n(u∗−a)n−1 n. The last condition is quadratic on aand we can rewrite it as: c1(d)a2+c2(d)a+c3(d)>0, where we have: c1(d) = d2m2, c2(d)=2dmu∗n(u∗)m n(u∗−a)n−1 n−d(m−1), c3(d)=(u∗)2d2(m−1)2−2dn(u∗)m n(u∗−a)n−1 n(m+ 1) + n2(u∗)2m n(u∗−a)2n−1 n. Thus, since we have a parabola opening upwards, the quadratic constraint gives two solutions: a<−c2(d)−pc2(d)2−4c1(d)c3(d) 2c1(d), or a>−c2(d) + pc2(d)2−4c1(d)c3(d) 2c1(d), and simplifying the expressions, we get: a<u∗ m (m−1) −n d(u∗)m n(u∗−a)n−1 n−2sn(u∗)m n(u∗−a)n−1 n d , (33) a>u∗ m (m−1) −n d(u∗)m n(u∗−a)n−1 n+ 2sn(u∗)m n(u∗−a)n−1 n d . (34) 46 Note that since d>1, conditions 4.3.1 and 4.3.1 do not contradict. Moreover, the constraint in (33) makes the condition 4.3.1 superfluous. However, (34) contradicts with condition 4.3.1, and cannot be satisfied. Therefore, we have two boundary curves (lower and upper) on the (b,a) plane which bound the Turing space. a>u∗ m(m−1) −n(u∗)m n(u∗−a)n−1 n, a<u∗ m (m−1) −n d(u∗)m n(u∗−a)n−1 n−2sn(u∗)m n(u∗−a)n−1 n d . It is easy to plot the above expression when n= 1 and m∈ {2, 3, 4}. Using u∗as a parameter going from 0 to ∞and using b=u∗−a. alower >u∗ m((m−1) −(u∗)m) , (35) aupper <u∗ m (m−1) −(u∗)m d−2r(u∗)m d!. (36) However, when n∈ {2, 3, 4}, rather than using the analytic expressions, we can numerically evaluate the code used to verify the Turing necessary conditions across a range of values. The resulting pairs (b,a) that satisfy the conditions are plotted, see Figure 18. The lower curve is independent of d, thus, the size of the Turing space can be changed by tuning d. The size of the Turing space is proportional to d, increasing dleads to a greater Turing space. Figures 17 and 18 show that b>ain most cases, as the Turing space mainly lies below the line a=b. Moreover, while the shape of the Turing space remains similar, its size increases for higher pairs of stoichiometric numbers, see Figure 18. Figure 17: Turing space for parameters (b,a) of the generalized Schnakenberg model (31) for different values of diffusion ratio d. The bounding curves are given by (35), (36). 47 Figure 18: Turing space for parameters (b,a) of the generalized Schnakenberg model (31) for d= 40. Simulations The implementation of System (31) follows the approach described in Section 3.2. Zero-flux boundary conditions are imposed to ensure no flux across the domain boundaries. The initial conditions for the receptor (u) and ligand (v) concentrations are set as small amplitude random perturbations around the homogeneous steady state (u∗,v∗) given in (32). Specifically, the initial conditions for Figures 19,20 and 21 are given by: u0= (a+b) + 0.01 ·rand(Nx), v0=b um 01 n + 0.01 ·rand(Nx). where Nxis the number of spatial points, aand bare non-negative constants and rand(Nx) generates uniformly distributed random values. These simulations aim to show if the different stoichiometric parameters mand ncan lead to the emergence of Turing patterns. We use arbitrary values that satisfy the Turing necessary conditions 4.3.1 for a,b(given in the figure captions below) and we set d= 40 and γ= 100. Figures 19,20 and 21 show that patterns are formed for all combinations of m∈ {2, 3, 4}and n∈ {1, 2, 3, 4}. Moreover, the peaks in receptor and ligand concentrations occur alternately, with the receptor concentration peaks being significantly larger and sharper than those of the ligand. 48 Figure 19: Generalized Schnakenberg model (31) for m= 2, n= 1, 2, 3, 4 and a= 0.125, b= 0.42 Figure 20: Generalized Schnakenberg model (31) for m= 3, n= 1, 2, 3, 4 and a= 0.2, b= 1.3. 49 Figure 21: Generalized Schnakenberg model (31) for m= 4, n= 1, 2, 3, 4 and a= 0.4, b= 1. 50 A. Codes 1import numpy as np 2import matplotlib . pyplot as plt 3 4# Parameters 5K2 = 0.012 #\ mu = 12 nM 6EP2_tot = 0.004 # \ muM = 4nM 7 8def EP2_star (x, n): 9return (x**n * EP2_tot ) / (x**n + K2) 10 11 12 x = np . linspace (0 , 1.5 , 1000) 13 14 plt . figure ( figsize =(8 , 5) ) 15 for nin [1, 2, 3, 4, 5]: 16 plt . plot (x, EP2_star (x , n) , label =f’n␣=␣{n}’) 17 18 # Labels and title 19 plt . title (r"$[\ text {EP2 }^*] $␣for ␣ different ␣ Hill ’s␣ coefficients ") 20 plt . xlabel (r" $PGE2$␣($\mu␣M$)") 21 plt . ylabel (r"$[\ text {EP2 }^*] $␣($\mu␣M$)") 22 plt . legend ( title =" Hill ’s␣ coefficient ") 23 plt. show () Listing 1: Python Script for Figure 4. 1import numpy as np 2import matplotlib . pyplot as plt 3from scipy . integrate import odeint 4 5# Parameters 6k1 = 5 #s^ -1 7k2 = 0.07 #s^-1 8k3 = 0.7 #\ muM ^-1 s^ -1 = 0.7 e+6 M^-1 s^-1 9k4 = 18.9e -3 #s^-1 10 beta_gamma_tot = 0.005 #\ muM = 5nM 11 alpha_s_tot = 2.3 #\muM 12 EP2_tot = 0.004 #\mu = 4nM 13 K1 = 0.8 #\ muM 14 K2 = 0.012 #\ mu = 12 nM 15 16 n=4# Hill ’s coefficient 17 18 k9 = 6.713 # s^ -1 19 k11 = 0.72 # s^-1 20 k12 = 6.85 # s^-1 21 K4 = 0.2 # \ muM = 200 nM 22 K6 = 0.09 # \ muM = 90 nM 23 K7 = 2.6 # \muM 24 K8 = 0.15 # \muM 57 25 AC_tot = 0.029 # \ muM = 28 nM 26 PDE4_tot = 0.115 # \muM = 115 nM 27 PDE3_tot = 0.0025 # \muM = 2.5 nM 28 C1 = 11 29 30 kw = 6.713 #s^-1 31 Kas = 0.2 #\muM 32 dw = 8.66 #s^-1 33 Kw = 1.21 #\ muM 34 A_tot = 0.0497 #\ muM 35 P_tot = 0.039 #\ muM 36 37 kd = 0.0058 # s ^-1 ( dissociation constant ) 38 ka = kd / K2 # Association constant 39 40 # Calculate EP2_star for each PGE2 value 41 PGE2_vals = [0.01 , 0.1 , 1, 10] 42 # PGE2_vals = np. linspace (0.1 ,1 ,10) 43 EP2_star_vals = [( PGE2 ** n * EP2_tot ) / ( PGE2 ** n + K2) for PGE2 in PGE2_vals ] 44 45 def ODE(vars, t): 46 beta_gamma , alpha_s_star , cAMP = vars 47 d_beta_gamma = (( k1/K1) * EP2_star + k4) * beta_gamma_tot - (( k1/K1 ) * EP2_star + k4 + k3 * alpha_s_tot - k3 * beta_gamma_tot ) * beta_gamma + k3 * beta_gamma * alpha_s_star - k3 * beta_gamma **2 48 d_alpha_s_star = (k1/K1) * EP2_star * beta_gamma_tot - (k1/K1) * EP2_star * beta_gamma - k2 * alpha_s_star 49 50 # Reference 1 51 # d_cAMP = k9 * (( AC_tot * alpha_s_star ) * (K6 + C1 * beta_gamma )) / (( alpha_s_star + K4) * (K6 + beta_gamma )) - (k11 * PDE4_tot * cAMP) / ( cAMP + K7 ) - (k12 * PDE3_tot * cAMP ) / ( cAMP + K8) 52 53 # simplified ( superactivation term out ) 54 # d_cAMP = k9 * ( AC_tot * alpha_s_star ) / ( alpha_s_star + K4) - (k11 * PDE4_tot * cAMP) / ( cAMP + K7) - (k12 * PDE3_tot * cAMP) / ( cAMP + K8 ) 55 56 # more simplified ( linear term ) 57 # d_cAMP = k9 * ( AC_tot * alpha_s_star ) / K4 - (k11 * PDE4_tot * cAMP ) / (K7 ) - (k12 * PDE3_tot * cAMP ) / (K8) 58 59 # Reference 2 60 d_cAMP = (kw * A_tot * alpha_s_star ) / ( alpha_s_star + Kas) - (dw * P_tot * cAMP ) / ( cAMP + Kw) 61 62 return [ d_beta_gamma , d_alpha_s_star , d_cAMP ] 63 64 # Time range for simulation 65 t = np . linspace (0 , 200 , 300) 66 58 67 # Initial conditions 68 ic = [6e -5, 0, 0] 69 70 # Titles and labels for the plots 71 plots = [ 72 {"i":0," title ":r"Trajectories␣of␣$[\ beta \ gamma ]$","ylabel":r"$[\ beta \ gamma ]$␣$(\ mu␣\ text {M}) $"}, 73 {"i":1," title ":r"Trajectories␣of␣$[\ alpha_s ^*] $","ylabel":r"$[\ alpha_s ^*]$␣$(\ mu ␣\ text {M }) $"}, 74 {"i":2," title ":r" Trajectories ␣of␣[ cAMP ]" ,"ylabel":r"[ cAMP]␣$(\ mu ␣\ text{M })$"} 75 ] 76 77 for iin plots : 78 plt . figure ( figsize =(10 , 6) ) 79 for PGE2 , EP2_star in zip ( PGE2_vals , EP2_star_vals ): 80 sol = odeint (ODE , ic , t) 81 plt . plot (t , sol [:, i["i"]] , label =f ’PGE2␣=␣{ PGE2 :.1f}’) 82 plt . title (i[" title "], fontsize =16) 83 plt . xlabel (" Time ␣(s)", fontsize =14) 84 plt . ylabel (i["ylabel"], fontsize =14) 85 #plt . grid (True) 86 plt . legend () 87 plt. show () 88 89 def ligand_dyn ( vars, t): 90 beta_gamma , alpha_s_star , PGE2 , EP2_star , EP2 , cAMP = vars 91 d_beta_gamma = (( k1/K1) * EP2_star + k4) * beta_gamma_tot - (( k1/K1 ) * EP2_star + k4 + k3 * alpha_s_tot - k3 * beta_gamma_tot ) * beta_gamma + k3 * beta_gamma * alpha_s_star - k3 * beta_gamma **2 92 93 d_alpha_s_star = (k1/K1) * EP2_star * beta_gamma_tot - (k1/K1) * EP2_star * beta_gamma - k2 * alpha_s_star 94 95 d_cAMP = (kw * alpha_s_star * A_tot ) / ( alpha_s_star + Kas ) - (dw * P_tot * cAMP ) / ( cAMP + Kw) 96 97 # Equations for ligand - receptor dynamics 98 d_PGE2 = kd * EP2_star - ka * PGE2 **n * EP2 99 d_EP2_star = ka * PGE2 **n * EP2 - kd * EP2_star 100 d_EP2 = kd * EP2_star - ka * PGE2 **n * EP2 101 return [ d_beta_gamma , d_alpha_s_star , d_PGE2 , d_EP2_star , d_EP2 , d_cAMP ] 102 103 # Time range for the simulation 104 t = np . linspace (0 , 200 , 300) 105 106 # Initial PGE2 values 107 PGE2_iv = [0.01 , 0.1 , 1, 10] 108 109 # Variables to plot and their corresponding indices in the solution array 110 plots = [ 111 {"i":0," title ":r"Trajectories␣of␣$[\ beta \ gamma ]$","ylabel":r"$[\ 59 beta \ gamma ]$␣$(\ mu␣M)$"}, 112 {"i":1," title ":r"Trajectories␣of␣$[\ alpha_s ^*] $","ylabel":r"$[\ alpha_s ^*]$␣$(\ mu ␣M)$"}, 113 {"i":2," title ":r" Trajectories ␣of␣[ PGE2 ]" ,"ylabel":r"[ PGE2]␣$(\ mu ␣M) $"}, 114 {"i":3," title ":r"Trajectories␣of␣$[\ text{ EP2 }^*]$","ylabel":r"$[\ text{ EP2 }^*] $␣$(\ mu ␣M)$"}, 115 {"i":4," title ":r" Trajectories ␣of␣[ EP2]" ,"ylabel":r"[EP2]␣$(\ mu ␣M)$" }, 116 {"i":5," title ":r" Trajectories ␣of␣[ cAMP ]" ,"ylabel":r"[ cAMP]␣$(\ mu ␣M) $"}, 117 ] 118 119 # Loop over each variable to generate plots 120 for iin plots : 121 plt . figure ( figsize =(10 , 6) ) 122 for PGE2_init in PGE2_iv: 123 # Set the initial conditions 124 ic = [6e -5, 0, PGE2_init , 0, EP2_tot , 0] 125 126 # Solve the system of ODEs 127 solution = odeint ( ligand_dyn , ic , t) 128 129 # Plot the solution for the current variable 130 plt. plot (t, solution [:, i["i"]] , label =f ’PGE2 ␣=␣{ PGE2_init :.2 f}’) 131 132 # Add labels , title , and legend 133 plt . title (i[" title "], fontsize =16) 134 plt . xlabel (" Time ␣(s)", fontsize =14) 135 plt . ylabel (i["ylabel"], fontsize =14) 136 plt . legend () 137 plt. show () Listing 2: Python Script for Figures 5and 6. 1import numpy as np 2from scipy . optimize import fsolve 3import matplotlib . pyplot as plt 4from scipy . sparse import diags 5from scipy.sparse.linalg import spsolve 6 7# ODE model of the ligand - receptor dynamics for the boundary fluxes of the PDE 8 9def backward_euler (PGE2 , EP2 , EP2_star , k_d , k_a , n, dt ): 10 def ligand_dyn ( vars): 11 PGE2_new , EP2_new , EP2_star_new = vars 12 U1 = k_d * EP2_star_new 13 U2 = k_a * ( PGE2_new **n) * EP2_new 14 return [ 15 PGE2_new - PGE2 - dt * (U1 - U2), 16 EP2_star_new - EP2_star - dt * (U2 - U1), 17 EP2_new - EP2 - dt * (U1 - U2) 60 18 ] 19 20 initial_condition = [ PGE2 , EP2 , EP2_star ] 21 solution = fsolve ( ligand_dyn , initial_condition ) 22 return solution 23 24 # Parameters 25 K2 = 0.012 #\ mu = 12 nM 26 k_d = 0.0058 # s^ -1 ( dissociation constant ) 27 k_a = k_d / K2 # s^ -1 ( association constant ) 28 n=4 # Hill coefficient 29 30 dt = 0.01 # Time step size 31 T = 200 # Total simulation time 32 33 num_steps = int(T / dt) 34 time_points = np. linspace (0, T, num_steps + 1) 35 PGE2_solutions = np . zeros ( num_steps + 1) 36 EP2_solutions = np . zeros ( num_steps + 1) 37 EP2_star_solutions = np . zeros ( num_steps + 1) 38 39 # Initial conditions 40 PGE2_0 = 10 # 0.01 ,0.1 ,1 ,10 41 PGE2_solutions[0] = PGE2_0 42 EP2_solutions [0] = 0.004 43 EP2_star_solutions [0] = 0 44 45 # Computing solutions 46 for step in range (0 , num_steps ): 47 PGE2_prev = PGE2_solutions [step ] 48 EP2_prev = EP2_solutions [ step] 49 EP2_star_prev = EP2_star_solutions [ step] 50 PGE2_new , EP2_new , EP2_star_new = backward_euler ( 51 PGE2_prev , EP2_prev , EP2_star_prev , k_d , k_a , n, dt 52 ) 53 PGE2_solutions [ step + 1] = PGE2_new 54 EP2_solutions [ step + 1] = EP2_new 55 EP2_star_solutions [ step + 1] = EP2_star_new 56 57 # Plotting the solutions 58 plt . figure ( figsize =(10 , 6) ) 59 plt . plot ( time_points , EP2_solutions , label ="EP2 ", linestyle ="-", linewidth =2) 60 plt . plot ( time_points , EP2_star_solutions , label =" EP2*", linestyle ="--", linewidth =2) 61 plt . xlabel (" Time ␣(s)") 62 plt . ylabel (r" Concentration ␣($\ mu$M)") 63 plt . title (" Trajectories ␣ of␣[ EP2 ]␣and ␣[ EP2 *]") 64 plt . legend () 65 plt. grid () 66 plt. show () 67 61 68 desired_time = 10.0 69 time_index = int(desired_time / dt) # Compute the time step index 70 EP2_star_value = EP2_star_solutions [ time_index ] 71 72 # Reaction - diffusion model 73 74 # Parameters 75 D1 = 1 76 D2 = 0.2 77 78 k1 = 5 79 k2 = 0.07 80 k3 = 0.7 81 k4 = 18.9e -3 82 K1 = 0.8 83 beta_gamma_tot = 0.005 84 alpha_s_tot = 2.3 85 EP2_tot = 0.004 86 n=4 87 88 kw = 6.713 89 Kas = 0.2 90 dw = 8.66 91 Kw = 1.21 92 A_tot = 0.0497 93 P_tot = 0.039 94 95 # Discretization Parameters 96 a=0 97 b = 70 98 T = 100 # Final time 99 dx = 0.01 100 dt = 0.01 101 Nx = int ((b - a) / dx) # Number of discretization points in space 102 Nt = int (T / dt) # Number of discretization points in time 103 x = np . linspace (a , b, Nx ) 104 t_lin = np. arange (0 , T + dt , dt) 105 106 # Define the matrix A obtained from discretization 107 108 lambd_1 = D1 * (dt / dx **2) 109 lambd_2 = D2 * (dt / dx **2) 110 111 # Helper function to construct matrices with zero -flux boundary conditions 112 def construct_matrix (Nx , lambd ): 113 upper = -lambd * np. ones(Nx -1) 114 main = (1 + 2 * lambd ) * np. ones(Nx) 115 lower = -lambd * np. ones(Nx -1) 116 diagonals = [lower , main , upper ] 117 A = diags ( diagonals , offsets =[ -1 , 0, 1] , format=" csr") 118 A[0, 0] = 1 + lambd # Zero - flux boundary condition at left 119 A[Nx -1 , Nx -1] = 1 + lambd # Zero - flux boundary condition at right 62 120 return A 121 122 # Construct matrices 123 A_beta_gamma = construct_matrix (Nx , lambd_2 ) 124 A_alpha_star = construct_matrix (Nx , lambd_2 ) 125 A_cAMP = construct_matrix (Nx , lambd_1 ) 126 127 # Solve the system of equations using ’spsolve ’ 128 129 # Initial conditions 130 beta_gamma_in = np .full (Nx , 6e -5) 131 alpha_star_in = np. full (Nx , 0.0) 132 cAMP_in = np. full(Nx , 0.0) 133 134 # Lists to store solution 135 store_beta_gamma = [] 136 store_alpha_star = [] 137 store_cAMP = [] 138 139 # Append initial conditions to the lists 140 store_beta_gamma . append ( beta_gamma_in . copy ()) 141 store_alpha_star . append ( alpha_star_in . copy ()) 142 store_cAMP . append ( cAMP_in . copy ()) 143 144 # Time stepping 145 for time_step_index , time_step in enumerate ( t_lin [: -1]) : # Exclude the last time point 146 EP2_star_next = EP2_star_solutions [ time_step_index + 1] 147 148 b_beta_gamma = beta_gamma_in + dt * ( 149 -k3 * beta_gamma_in * ( beta_gamma_in - alpha_star_in - beta_gamma_tot + alpha_s_tot ) 150 + k4 * (beta_gamma_tot - beta_gamma_in)) 151 # Update boundaries in b vector for beta_gamma 152 b_beta_gamma [0] += (dt / dx) * (k1 / K1) * EP2_star_next * ( beta_gamma_tot - beta_gamma_in[0]) 153 b_beta_gamma [ -1] += (dt / dx) * (k1 / K1) * EP2_star_next * ( beta_gamma_tot - beta_gamma_in[-1]) 154 155 b_alpha_star = alpha_star_in + dt * (-k2 * alpha_star_in ) 156 # Update boundaries in b vector for alpha_star 157 b_alpha_star [0] += (dt / dx) * (k1 / K1) * EP2_star_next * ( beta_gamma_tot - beta_gamma_in[0]) 158 b_alpha_star [ -1] += (dt / dx) * (k1 / K1) * EP2_star_next * ( beta_gamma_tot - beta_gamma_in[-1]) 159 160 b_cAMP = cAMP_in + dt * (-dw * ( P_tot * cAMP_in ) / ( cAMP_in + Kw)) 161 # Update boundaries in b vector for cAMP 162 b_cAMP [0] += kw * ( A_tot * alpha_star_in [0]) / ( alpha_star_in [0] + Kas ) 163 b_cAMP [ -1] += kw * ( A_tot * alpha_star_in [ -1]) / ( alpha_star_in [ -1] + Kas) 164 63 165 # Solve the linear system Au = b for each variable 166 beta_gamma = spsolve ( A_beta_gamma , b_beta_gamma ) 167 alpha_star = spsolve ( A_alpha_star , b_alpha_star ) 168 cAMP = spsolve ( A_cAMP , b_cAMP ) 169 170 # Update solutions 171 beta_gamma_in = beta_gamma 172 alpha_star_in = alpha_star 173 cAMP_in = cAMP 174 175 # Update solution 176 store_beta_gamma . append ( beta_gamma . copy ()) 177 store_alpha_star . append ( alpha_star . copy ()) 178 store_cAMP . append ( cAMP . copy ()) 179 180 # Plot solutions over space 181 182 fig , ax = plt. subplots (1 , 2, figsize =(22 , 8) ) 183 184 # Plot the initial conditions ( time step 0) 185 ax [0]. plot (x , store_beta_gamma [0] , ’b-’, label =r’$[\ beta \ gamma ] _0$’, linewidth =3) 186 ax [0]. plot (x , store_alpha_star [0] , ’m-’, label =r’$[\ alpha_s ^*] _0$’, linewidth =3) 187 ax [0]. plot (x , store_cAMP [0] , ’y-’, label =r’[cAMP ] $_0$’, linewidth =3) 188 ax [0]. set_xlabel (’x’, fontsize =25) 189 ax [0]. set_title (" Initial ␣ Conditions ", fontsize =20) 190 191 # Plot the final conditions (last time step) 192 ax [1]. plot(x, store_beta_gamma [ -1], ’b-’, label =fr’$[\ beta \ gamma ]$(T ={T }) ’, linewidth =3) 193 ax [1]. plot(x, store_alpha_star [ -1], ’m-’, label =fr’$[\ alpha_s ^*] $(T ={ T}) ’, linewidth =3) 194 ax [1]. plot(x, store_cAMP [ -1] , ’y-’, label =fr ’[ cAMP ]( T={T }) ’, linewidth =3) 195 ax [1]. set_xlabel (’x’, fontsize =25) 196 ax [1]. set_title (f" Solution ␣at ␣T ={T}" , fontsize =20) 197 198 fig . suptitle (fr" Reaction - Diffusion ␣ System ␣for ␣[PGE2 ] $_0$={ PGE2_0 }␣$\ mu$M", fontsize =20) 199 for iin range (2) : 200 plt. sca (ax[i]) 201 plt . xticks ( fontsize =20 , family =’serif ’) 202 plt . yticks ( fontsize =20 , family =’serif ’) 203 ax [i ]. tick_params ( axis =’both ’, which = ’major ’, length =8) 204 ax [i ]. tick_params ( axis =’both ’, which = ’minor ’, length =4) 205 ax [i]. legend ( loc =’best ’, fontsize =20) 206 207 plt. show () Listing 3: Python Script for Figures 8,9and 10. 1import numpy as np 2from scipy . sparse import diags 64 3from scipy.sparse.linalg import spsolve 4import matplotlib . pyplot as plt 5 6# Reaction - diffusion model 7# Parameters 8D1 = 1 9D2 = 0.2 10 D3 = 0.08 # 0.00208 #0.2 PONER 0.05!!! 11 12 k1 = 5 13 k2 = 0.07 14 k3 = 0.7 15 k4 = 18.9e -3 16 K1 = 0.8 17 beta_gamma_tot = 0.005 18 alpha_s_tot = 2.3 19 # EP2_tot = 0.004 20 n=4 21 22 K2 = 0.012 23 kd = 0.0058 24 ka = kd / K2 25 26 kw = 6.713 27 Kas = 0.2 28 dw = 8.66 29 Kw = 1.21 30 A_tot = 0.0497 31 P_tot = 0.039 32 33 # Discretization Parameters 34 a=0 35 b = 70 36 T = 100 # Final time 37 dx = 0.01 38 dt = 0.01 39 Nx = int ((b - a) / dx) # Number of discretization points in space 40 Nt = int (T/dt) # Number of discretization points in time 41 x = np . linspace (a , b, Nx ) 42 t_lin = np. arange (0 , T + dt , dt) 43 44 # Define the A matrix obtained from discretization 45 46 lambd_1 = D1 * (dt / dx **2) 47 lambd_2 = D2 * (dt / dx **2) 48 lambd_3 = D3 * (dt / dx **2) 49 50 # Helper function to construct matrices with zero -flux boundary conditions 51 def construct_matrix (Nx , lambd ): 52 upper = -lambd * np. ones(Nx -1) 53 main = (1 + 2 * lambd ) * np. ones(Nx) 54 lower = -lambd * np. ones(Nx -1) 65 55 diagonals = [lower , main , upper ] 56 A = diags ( diagonals , offsets =[ -1 , 0, 1] , format=" csr") 57 A[0, 0] = 1 + lambd # Zero - flux boundary condition at left 58 A[Nx -1 , Nx -1] = 1 + lambd # Zero - flux boundary condition at right 59 return A 60 61 # Define matrices for each variable 62 A_1 = construct_matrix (Nx , lambd_1 ) # PGE2 63 A_2 = construct_matrix (Nx , lambd_3 ) # EP2 * 64 A_3 = construct_matrix (Nx , lambd_3 ) # EP2 65 A_4 = construct_matrix (Nx , lambd_2 ) # beta_gamma 66 A_5 = construct_matrix (Nx , lambd_2 ) # alpha_s_star 67 A_6 = construct_matrix (Nx , lambd_1 ) # cAMP 68 69 # Initial condition for EP2_star 70 def ic_EP2 (x , Nx , EP2_tot , num_clusters , sigma ): 71 EP2_in = np. zeros (Nx ) # Initialize with zeros 72 cluster_centers = np . linspace (x [0] , x[ -1] , num_clusters + 2) [1: -1] 73 74 for center in cluster_centers : 75 EP2_in += EP2_tot * np. exp (-(x - center)**2 / (2 * sigma **2) ) # Set peak height to EP2_tot 76 77 return EP2_in 78 79 EP2_tot = 0.004 # Peak height equal to EP2_tot 80 num_clusters = 2 # Number of clusters 81 sigma = 5.0 # Width 82 83 EP2_in = ic_EP2 (x, Nx , EP2_tot , num_clusters , sigma) 84 85 # Solve the system of equations using ’spsolve ’ 86 # Initial conditions 87 PGE2_0 = 0.1 88 PGE2_in = np. full (Nx , PGE2_0 ) 89 EP2_star_in = np. full (Nx , 0.0) 90 EP2_tot = 0.004 91 EP2_in = np. full (Nx , EP2_tot ) 92 beta_gamma_in = np .full (Nx , 6e -5) 93 alpha_star_in = np. full (Nx , 0.0) 94 cAMP_in = np. full(Nx , 0.0) 95 96 # Lists to store solution 97 store_PGE2 = [] 98 store_EP2_star = [] 99 store_EP2 = [] 100 store_beta_gamma = [] 101 store_alpha_star = [] 102 store_cAMP = [] 103 104 # Append initial conditions to the lists 105 store_PGE2 . append ( PGE2_in ) 66 79 plt . xticks ( fontsize =20 , family =’serif ’) 80 plt . yticks ( fontsize =20 , family =’serif ’) 81 ax [i ]. tick_params ( axis =’both ’, which = ’major ’, length =8) 82 ax [i ]. tick_params ( axis =’both ’, which = ’minor ’, length =4) 83 ax [i]. legend ( loc =’best ’, fontsize =20) 84 85 plt. show () Listing 8: Python Script for Figures 19,20 and 21. 73