MATHEMATICAL ANALYSIS AND NUMERICAL SIMULATION OF A MODEL OF MORPHOGENESIS

We consider a simple mathematical model of distribution of morphogens (signaling molecules responsible for the differentiation of cells and the creation of tissue patterns). The mathematical model is a particular case of the model proposed by Lander, Nie and Wan in 2006 and similar to the model presented in Lander, Nie, Vargas and Wan 2005. The model consists of a system of three equations: a PDE of parabolic type with dynamical boundary conditions modelling the distribution of free morphogens and two ODEs describing the evolution of bound and free receptors. Three biological processes are taken into account: diffusion, degradation and reversible binding. We study the stationary solutions and the evolution problem. Numerical simulations show the behavior of the solution depending on the values of the parameters.


1.
Introduction.M orphogenesis (the creation "genesis" of shapes "morphe") has been studied from the early 20th century, but only in recent years, growth factors have been identified as morphogens.The formation of the embryo can be understood only as a global process and therefore global phenomena as differentiation of tissues, formation of organs and its organization have to be considered to understand it.The differentiation of the cell, the key process on the formation of embryos, depends on its position.The cell receives the information of its position by measuring the concentration of signaling molecules, named morphogens.
Morphogens are synthesized at signaling localized sites and spread into the body creating gradients in its concentration, otherwise, a constant distribution of morphogens would create an homogeneous differentiation of cells.How the gradients arise is an unclear and controversial question and central issue in Development Biology.Theoretical and experimental scientists consider two main theories to explain the formation of gradients of morphogens: diffusion theory, where morphogens are spread by diffusion through the extracellular matrix and the positional theory (see Kerszberg and Wolpert [7]) which suggests that morphogen is transferred by contact.
Once the morphogens arrive to the cell surface they bind to receptors and other kind of molecules.The diffusion theory considers slow degradation of products and reversible binding (see Lander, Nie, Wan [10]) in contrast the positional theory does not consider degradation (see Kerszberg and Wolpert [7]).
Lander, Nie and Wan [10] studied numerically several mathematical models and focused on the Drosophila wing disc.They obtain that diffusion process and nonhomogeneous distribution of morphogen are not in contradiction.Tello [18] studied the mathematical model and the asymptotic behavior of the solutions of the model proposed in Lander, Nie and Wan [10].
Lander, Nie, Vargas and Wan [9] and Lander, Nie and Wan [11] proposed several models of differential equations.The models consider a PDE of parabolic type to describe the evolution of morphogens and several ODEs to model the receptor and the bound-receptor.They study the steady states and the linear stability of them under the action of a source in a region of the domain.In the present work, we shall analyze the evolution problems, corresponding to some particular cases of the model proposed by Lander, Nie, Vargas and Wan [9].
Merkin and Sleeman [15] have studied the system proposed by Lander, Nie and Wan [10] with degradation and without it.They provide an analysis of the models under the assumption of constant concentration of morphogens at the boundary x = 0 and gradient of morphogens equals to zero at infinity.The authors prove that the case where the bound morphogen complex is not degraded, the free morphogen profile is essentially linear and spreads as a square root law.
Recently Merkin, Needham and Sleeman [14] have introduced a chemosensitivity term in the model to describe morphogen concentration.They have presented results on the existence and uniqueness of classical solutions and self-similarity.Their numerical simulations have showed periodic pulse solutions.
Lou, Nie and Wan [13] consider a model with two species of morphogens, Dpp and Sog.The system consists of three PDEs of parabolic type and one ODE.The authors study the steady states of the system and prove the existence of nontrivial gradients in the solution under biologically meaningful assumptions.Numerical simulations of the evolution of the system show the behavior of the solutions for different Sog production rates.For high Sog production rates, the solution exhibits an intense concentration of the complex Dpp-receptor in the region of Dpp production.The magnitude of the concentration increases as Sog production rate increases.
In this work we study the case of diffusive transport of morphogens.In Section 2 we consider the mathematical model proposed by Lander, Nie and Wan [11], the simplifications we shall consider and the resulting mathematical model to be treated throughout the present work.The particular case does not consider the effects of the processes in the interior of the cell.The resulting model is similar to the model described by Lander, Nie, Vargas and Wan [9].Section 4 is devoted to the steady states.The boundary conditions are studied in Section 5 and the results concerning existence and uniqueness of solutions are introduced in Section 6.Finally, Section 7 is about the numerical resolution of the model.We present some numerical results and simulations of the model for some particular values of the parameters and enumerate some conclusions which can be derived from the results.The description of the scheme of resolution is presented in the Appendix.
2. The mathematical model.Different models of distribution of morphogens have been introduced by several authors in the last decade.We study, in the following sections, a particular case of a mathematical model proposed by Lander, Nie and Wan [11].They considered a 1D reaction-diffusion system of partial differential equations and auxiliary conditions governing the morphogen activities along the anterior-posterior axis of a Drosophila wing disc.We consider below the mathematical model proposed by Lander, Nie and Wan [11], see also Lander, Nie, Vargas and Wan [9] for the details of the modeling.The unknowns of the problem represent the normalized concentrations of the following species: • u is the concentration of morphogen, • b is the concentration of ligand-receptor complexes bound to cell surface membrane, • c is the concentration of ligand-receptor complexes in the cell interior • d is the extracellular receptor concentration, • e is the concentration of receptors in the interior of the cell.
The system of equations is defined in I := (0, 1) and with the boundary condition x = 1 : u = 0 t > 0, (7) and the initial data where h 0 , f 0 , g 0 , j 0 , j 1 , k 0 , k 1 , ν, σ 0 , e 0 and ω 0 are parameters of the problem.The partial differential equation (1) governs the rate of change of morphogen concentration "u", h 0 and f 0 are the binding and disassociation rate constants.Equations ( 2) and (3) model the processes of endocytosis (process by which cells absorb material) and exocytosis (process by which cells secret components to the extracellular matrix) and consider degradation of the ligand-receptor complexes.
The evolution of the receptor concentrations, in the intracellular and in the extracellular, are considered in equations ( 4) and (5).Formation, dissociation and degradation of morphogen-receptor complexes, degradation of receptors, synthesis of new receptors are also considered in (4) and (5).
Lander et al. [11] focus their studies on the analysis of the steady state corresponding to this model, considering both cases, σ 0 = 0 and σ 0 > 0. For simplification we consider throughout the paper the particular case We assume that the influence of the concentration of the ligand-receptor complex in the interior of the cell is not relevant for the evolution of the complex in the cell surface membrane, i.e. k 0 = 0, (10) Then, equation ( 3) is uncoupled and c is obtained after integration To simplify, we consider the concentration of receptors in the interior of the cell is constant on time, i.e.
(11) The system is reduced to a system of three equations with dynamical boundary conditions.We introduce the new unknowns and the parameters Then, under assumptions ( 9)- (11) the system (1)-( 8) becomes with boundary conditions ∂u ∂t = ν − uw + v, at x = 0 for t > 0, u = 0 at x = 1 (15) and initial data u = v = 0, w = 1 at t = 0, x ∈ I .
In section 5 we study the boundary conditions.We will prove that the dynamical system at x = 0 presents a unique steady state for κ = ν, infinitely many for κ = ν = 0 and no steady states for κ = ν = 0.
In sections 4 and 5 we distinguish different cases depending on the parameters κ and ν: (17) represents the case where the production rate κ of receptors in the cell is larger than the production rate at x = 0 of the morphogen.( 18) and ( 19) consider the case where the balance is zero and negative respectively.
Assumptions (17) guarantees the non negativity of the values of the concentrations and therefore we exclude the non-biological admissible solutions.
3. Known results on functional analysis.In order to study the system of differential equations we introduce some functional spaces frequently used to solve partial differential equations.
Let Ω ⊂ IR n be an open and bounded set with regular boundary.Then we define the following spaces: ) are Hilbert Spaces .Let X(Ω) be a Banach space and • X be the norm defined in X(Ω) then The proof of Theorem 6.1 is based in Schauder fixed point theorem.In order to obtain the needed estimates we used several inequalities, listed below.
• Young s Inequality.Let a and b be nonnegative real numbers and p and p positive real numbers satisfying 1 ≤ p ≤ p ≤ ∞ with 1/p + 1/p = 1 then • Hölder s Inequality.Let p and p , positive numbers satisfying 1 ≤ p ≤ p ≤ ∞ with 1/p + 1/p = 1 and let f and g be real functions defined over Ω, then Notice that for Ω = I = (0, 1) we have that See for instance Gilbarg and Trudinger [5] for more details and proofs.
4. The steady states.The steady states of the problem are given by the system: with the boundary conditions: We combine (21) with ( 22) to obtain then equations ( 20)-( 22) become We introduce the normalized unknown ũ for simplicity, we drop the tilde and the equation (24) becomes (26) has been also obtained in Lou, Nie and Wan [13], for the time-independent solutions of a system considering blind-and free-morphogens.Explicit solutions may be found with an unknown parameter (c 1 ) implicitly defined (see equation ( 31)) as suggested in [13] which gives a range of admissible boundary condition.
However it is necessary to analyze the integral (31) to obtain that the studied boundary conditions are admissible.Uniqueness of solutions is a byproduct of this analysis.
Lander, Nie, Vargas and Wan [9] have consider the problem with a source term and nonlinear mixed boundary conditions.The problem with nonlinear boundary condition is more complicated and a monotone method is used to prove the existence of solutions.In [9], the existence and uniqueness of solutions is obtained by using a monotonicity method of Amman and Sattinger based on upper and lower solutions.The proof presented below is slightly more simpler and self-contained.We only use basic integration formulas and limits.The monotonicity of solutions is also studied in [13] and [9].
If κ = 0 and thanks to (23) we get hence, the boundary conditions for (26) for κ = 0 are the following • If κ = ν = 0 there is no solution to (21)-( 23) in x = 0. Notice that in the case κ = 0 where assumption (17) is satisfied, we have non negativity of the values of the concentrations.When assumption (17) is not satisfied we obtain non-biological admissible boundary conditions.
The case κ = 0 is simpler and two different options are possible: Since u(1) = 0, the solution is given by = 0 for any u(0), and the solutions are given by u(x) = A(1 − x), x ∈ I, for any positive constant A. The main result of this section is the following proposition: Proposition 3.For every k = ν there exists a unique solution "u" to (26), (27).Proof.We consider two different cases: Case I, α = ν κ−ν > 0. We multiply equation (26) by the negative part of u, then after integration by parts we obtain that any solution satisfies u ≥ 0 on I. Uniqueness of solutions is a consequence of the monotonicity and regularity of the function u u+1 for u ≥ 0.
The case α = 0 is the trivial case with a unique solution u = 0. We introduce the system of ODEs Then, after integration we obtain 1 2 and hence u (2ku − 2k ln(u + 1) + c 1 ) The case u has been excluded because there is no solution satisfying u > 0, u(0) = α > 0 and u(1) = 0. Notice that by mean value theorem, if such solution u exists, it has to satisfy u (x 0 ) = u(1) − u(0) = − α < 0 for some x 0 ∈ (0, 1), since u > 0, we have to exclude such case.Let F (c 1 , k, α) be defined by (31) then, after integration in (30), we have that there exists a solution if and only if there exists a constant c 1 such that Since F is a continuous and monotone function in c 1 , we have the existence of a unique c 1 > 0 for every α > 0 and therefore the existence of a unique solution u, implicitly defined by u α ds (2ks − 2k ln(s + 1) + c 1 ) Notice that from (30) we deduce that the steady state u is monotone decreasing.Notice also that if c 1 < 0, the solution satisfies u(1) > 0 which contradicts (6).Case II, α = ν κ−ν < 0. We multiply equation (26) by the positive part of u and after integration by parts we obtain that u ≤ 0 on I. Since u represents the concentration of morphogens, the solutions for κ − ν < 0 are non-biological admissible solutions.To complete the mathematical analysis of the problem we consider the solutions to (26) for α < 0.
The analysis is similar to case I, and u satisfies u (2ku − 2k ln |u + 1| + c 1 ) and F is defined by for c 1 > c * 1 = −2k α + 2k ln(− α + 1) and k > 0. The proof ends in the same way as previous case.2

5.
The boundary condition.We study in this section the boundary condition of the system.For reader s convenience we denote the variable (u(0, t), v(0, t), w(0, t) by (u, v, w).Then the boundary conditions at x = 0 are described by and the initial data Lemma 5.1.If κ = ν there exists a unique steady state of the system (33)-( 35), given by If κ = ν = 0 there is no steady states.If κ = ν = 0 there exists infinitely many steady states given by (u, 0, 0).
Proof.If κ = ν the steady states of the problem are given by the solutions of the equations Since u w − v = ν, we have that ν − µv = 0 and κ − ν − ηw = 0 which implies and hence If κ = ν, we subtract equation (37) to (39) and we obtain w = 0. We substitute in (38) to obtain v = 0. Then equations (37)-(39) have solutions if and only if κ = ν = 0 which proves the case κ = ν. 2 Remark.Notice that the positivity of the steady state is a direct consequence of assumption (17).
In order to study the existence of solutions for t ∈ (0, ∞) we first obtain some a priori estimates.Lemma 5.2.u, v and w satisfy Proof.The result is a consequence of the following facts Lemma 5.3.v and w are uniformly bounded.
Proof.We consider equations (34), (35) and the functions v f0 and w h0 , such that Then we have after integration that , which implies that v and w are uniformly bounded.
Proof.We divide equation (34) by f 0 and add to equation (33) to obtain After integration we get that which proves the Lemma.
In order to apply a fixed point argument of Schauder type, we introduce the functional space H 1 0 (I) and L 2 (0, T : H 1 0 (I)) defined in Section 3. Let T ∈ (0, ∞) and k defined by and A defined by Notice that since u(1, t) = 0, the norm L 2 (0, T : H 1 (I)) in A is equivalent to the norm in L 2 (0, T : H 1 0 (I)).We consider the functional J : A ⊂ L 2 (0, T : H 1 (I)) → L 2 (0, T : H 1 (I)) defined by J(ũ) = u where u is the solution to the problem with the boundary conditions: and the initial data: Before presenting the proof of the Theorem 6.1 we obtain the necessary a priori estimates.
Lemma 6.2.Let ũ ∈ A, then, the solutions v and w to (42), ( 43) and (45) satisfy Proof.(42), ( 43) is a linear system of two equations with bounded coefficients and therefore there exists a unique solution (v, w), moreover, it is bounded and for fixed x ∈ I the solutions v, w ∈ C 0 ([0, T ]).We introduce the functions ψ , ψ, Ψ and Ψ : IR → IR defined by Notice that Ψ = ψ and lim We multiply equation (42) by ψ (v) f0 , and (43) by ψ (w) h0 and integrate over (0, 1) to obtain: We take limits when → 0 to obtain Since ũ ≥ 0 and κ ≥ 0 we have that wũψ(v) ≤ ũΨ(w); vψ(w) ≤ Ψ(v); κψ(w) ≤ 0, and then We add above expressions to get We divide equation (42) by f 0 and equation ( 43) by h 0 , adding both expressions, we have we have, after integration that , (48) and ( 49) prove the Lemma. 2 Lemma 6.3.Let ũ ∈ A and v, w the solutions to (42), ( 43) satisfying (45) then, the solution u to (41), ( 44) with the initial data given by (45) satisfies Proof.Since u satisfies the equation ( 41), by Lemma 6.2, we have that we multiply (50) by (u−k t) + (where (•) + is the positive part function) and integrate by parts to obtain we have (u − k t) + = 0 at x = 0, 1 and therefore By Gronwall s lemma we have u ≤ k t.In the same way we multiply equation (41) by ψ (u) integrate by parts to obtain d dt Again, by applying Gronwall s lemma we have Ψ(u) = 0 which ends the proof of the Lemma. 2 Proof of Theorem 6.1.In order to obtain the existence of a fixed point of J we consider the following: I.-J is well defined.Let T < ∞.Since u ∈ C ∞ (0, T ) and w, v ∈ L ∞ (I × (0, T )) we have that there exists a unique solution u to (41) satisfying II.-J is continuous.Let ũi ∈ A for i = 1, 2, u i := J(ũ i ), v i and w i the solution to (42) and ( 43) for ũ = ũi .Then We multiply (54) by (u 1 − u 2 ) and integrate over I to obtain 1 2 applying Holder and Young inequalities and Lemma 6.3, it results 1 2 We take After integration in (57) over (0, T ) it results in particular we deduce 1 2 We integrate (57) over (0, t) to obtain and ( 59) over (0, T ) to have we have We introduce the constant k 4 := k 3 (1 + T ) and we combine (58) and (60) to obtain We take (v 1 − v 2 ) and (w 1 − w 2 ) as test functions in (55) and ( 56) and apply Holder and Young inequalities to get for a positive constant k 5 .( 62) and (63) imply after integration over (0, t) we have We integrate (64) over (0, T ) to obtain (61) and (65) prove the continuity of J.
It is a consequence of (53) and Lemma 6.3.
Notice that B 0 ⊂ B ⊂ B 1 and B 0 ⊂ B is a compact imbedding and is a Banach space.Let u, v and w the solutions to (41)-(45) for ũ ∈ A, we define f (x, t) := −uw + v, then, by Lemma 6.2 and Lemma 6.3 we have that Lemma 5.5, Lemma 5.2 and (66) imply where k 7 is a continuous function of T and defined for T < ∞ (see for instance Friedman [4]).Since the imbedding is compact (see J.L. Lions [12] Theorem 5.1 p. 58) we have that J(A) is a precompact subset of A.
We now apply Schauder's fixed point Theorem to Ã and J defined by Ã := {u ∈ L 2 (0, T : Thanks to I-IV there exists a fixed point of J in Ã and consequently a fixed point of J in A which proves the existence of at least one solution to the problem ( 12) -( 16).( 61) and (65) prove the uniqueness of solutions for T small enough and therefore uniqueness of solutions of the problem for all T < ∞. 2 7. Numerical results.In this section we shall present some numerical results obtained with different parametric values in the numerical resolution of the problem ( 12) - (16).In doing this we intent to show how the parameters determine the qualitative behavior of the three variables.The description of the scheme of resolution of the problem ( 12) -( 16) is presented in the Appendix.The fact that we are dealing with evolution problems, suggests to consider an iterative decoupling scheme.We shall solve separately in each step (in time) of the iterative scheme, a stationary problem for each of the three variables, the concentration of morphogens "u", the concentration of ligand-receptor complexes "v" and the concentration of receptors "w".The results show the influence of the parameters in the evolution of the variables u, v and w, from the quantitative and qualitative point of view.Some of the conclusions that we deduce from the computations are commented in the captions below Figure 2, Figure 3 and Figure 4. We also solve numerically the initial value problem at x = 0 defined by ( 13) -( 16) for some parameters satisfying the constraints κ − ν = 0, κ − ν > 0 and κ − ν < 0 respectively.In order to do this, we use a finite difference iterative scheme of step in time dt = 10 −6 .To precise, we solve first the evolution problem for u, considering an explicit Euler scheme.Next we solve the problem for v using again an explicit Euler scheme where the value of u has been already update.Finally we use the same solver for the problem for w, considering the updated values of v and u.In Figure 1 we present the graphics obtained for the three variables for some parametric values.8. Conclusions.We consider a mathematical model of distribution of morphogens, where morphogens are secreted in signaling points and transported by diffusion in the extracellular matrix.Mathematical analysis shows that gradients of morphogens may be produced by diffusion.We study a particular case of a mathematical model proposed by Lander et al [9] and [10].To simplify the problem we consider that the concentration of free receptors at the interior of the cells is constant on time.The problem is reduced to a system of three equations with dynamical boundary conditions.The behavior of the system is given by the evolution of morphogens at the secretion point x = 0 which stabilizes when the rate of production of receptors at the surface of the cell is larger than the rate of production of morphogens at x = 0. Small perturbations in the initial data will not produce any change in the final result.The other case (rate production of morphogens larger than the production of receptors) gives unbounded solutions with no Biological meaning.Several mathematical questions remain open concerning the behavior of the boundary conditions: Global boundedness for κ > ν and stability of the steady states which have been shown numerically.
The existence of steady states is studied when the parameters κ ν are positive.The proof given here is much simpler than the one presented in Lander, Nie, Vargas and Wan [9].
When the rate of production of receptors at the surface of the cell is larger than the rate of production of morphogens at x = 0, i.e. under assumption κ − ν > 0, the existence of positive steady state is guaranteed, otherwise negative steady state may appear.
Numerical simulations show an intensive concentration of morphogens and complex morphogen-receptor near the signaling point x = 0.The slope of the gradient of morphogens is mainly determined by the evolution of the morphogen at x = 0 which fast goes close to a constant value ( ν κ−ν ) and the distribution of free receptors remains almost homogeneous in space from the beginning of the experiment (see Figures 2, 3 and 4).The consequence of this distribution is a clear differentiations of cells close to the signaling point different than the one in the other part of domain.The large slope of the distribution of complex gives small measure for the amount of cells at the threshold value.The numerical simulations show that, for κ > ν, the solution tends to the steady state as fast as the concentration of morphogens at the signaling point stabilizes in the steady state given in Lemma 5.1.The system at the signaling point x = 0 is also analyzed numerically, where the solution stabilizes in a constant value (for κ > ν) and goes to ∞ for κ = ν = 1.In the case κ < ν the solution u goes to ∞ (see Lemma 5.5).
As a final conclusion we want to point out the relevance of the balance between the production of morphogens and the production of receptors, which determines the behavior of the system.We have simplified the model in Section 2 in order to focus in the process in the extracellular matrix.It is clear that the production of receptors at the surfaces cell depends on more complex mechanisms in the interior of the cell than those considered here where the production rate ν is assumed to be constant.Next, we proceed to solve the problem (68) and (71) with respect to the spatial variable x.We employ a finite element method which we describe below.Notice that in order to solve (69) and (70), we consider, after discretizing the spatial domain, the following equalities at each node x i of the mesh: We introduce the definition of weak solution of the problem.
Notice that any classical solution of (68) and ( 71) is also a weak solution.
Discretization in space.We consider an uniform mesh in [0, 1] consisting of the nodes M = {x 1 , x 2 , ....x N }, 0 = x 1 < x 2 ... < x N −1 < x N = 1 of size h = 1/N − 1.We use linear finite elements L n , n = 1, ... N associated to the mesh M. The discretized version of the problem is the following: The steps in time and space are, respectively dt = 10 −6 and h = 10 −3 .The system is solved using the Gauss Seidel method.The number of iterations depends on a fixed tolerance regarding the relative error.In doing so, we guarantee the convergence of the iterative scheme.In each step, the iterative scheme is started with the final value obtained in the previous step in time.The nodal values of the other two variables w and v, are computed via the identities given in (72) and (73).

Figure 1 .
Figure 1.The variables u, v and w at x = 0 in the graphics are represented by blue crosses, red diamonds and green circles, respectively.On the left, we present the results obtained for the parameters κ = 1, ν = 1, µ = 1 and η = 1.On the right, the values considered are κ = 10, ν = 1, µ = 1 and η = 1.