Intracellular protein dynamics as a mathematical problem

In this paper we undertake a mathematical analysis of a model of intracellular protein dynamics, i.e. protein and mRNA transport inside a cell, proposed in [1]. The model takes into account di ﬀ usive transport in the nucleus and cytoplasm, as well as active transport of protein molecules along microtubules in the cytoplasm. The model reproduces, at least in numerical simulations, the oscillatory changes in protein concentration observed in the experimental data. To our knowledge this is the ﬁrst paper that, in the multidimensional case, deals with a rigorous mathematical analysis of a model of intracellular dynamics with active transport on microtubules. In particular, in the present paper we prove well–posedness of the model in any space dimension. The model is a complex system of nonlinear PDEs with speciﬁc boundary conditions. It may be adapted to other signaling pathways.


Intracellular protein dynamics as a mathematical problem
Mirosław Lachowicz a , Martin Parisot b , Zuzanna Szymańska c

Introduction
Mathematical modeling of intracellular proteins dynamics is one of the most commonly undertaken challenges in mathematical biology.Of course, this follows from our belief that increased knowledge about these processes will allow us to better understand the functioning of a cell, especially if it is distorted, as it is in the case of many diseases like cancer.In the literature one can find many models describing the intracellular dynamics of different proteins, e.g.[2,3,4,5,6,7].However, most of these models, either ignores the spatial aspect, or treat it in a very simplified way, taking into account only the diffusion as a transport mechanism of molecules.Only recently new models that take into account directed transport that is held in the cytoplasm were published [1,8,9,10].
One of few such models, was proposed in 2014 by Szymańska at al., the model of a simple intracellular protein concentration control system, that took into account active transport along microtubules [1].The work focused not on particular protein dynamics but on the proper description of intracellular transport processes.The main goal was to provide a description that is fuller in terms of biology and correct from the mathematical point of view.Therefore the model is generic, however, corresponds to protein Hes1 simple regulatory system.
In the present paper, we show a well-posedness of a model of intracellular transport processes generalizing the system presented in [1] to the case of N couples mRNA/protein.The mathematical analysis is perform for an arbitrary space dimension and with non-linear parameters satisfying classical mathematical properties.To our knowledge, this is the first work that in multidimensional case deals with rigorous mathematical analysis of the model of intracellular protein dynamics including active transport along microtubules.The model can be describe as a system of hyperbolic and parabolic equation coupled through non-linear parameters and boundary condition.The mathematical analysis of such system is not trivial and we propose a analysis based on classical results for, on the one hand linear hyperbolic equation, on the second linear parabolic equation.The key point of the proof of well-posedness for the coupled nonlinear system lies in a priori estimations leading to contraction argument.Last but not least, some properties such as global evolution and positivity are given to ensure the biological relevance of the model.
For other paper with mathematical analysis of a gene regulatory network model in a 1-dimensional domain including linear stability analysis of the steady states we refer the reader to recent work of Chaplain at al. [11].The authors showed that a diffusion coefficient acts as a bifurcation parameter and give rise to a Hopf bifurcation.They however did not study directly the influence of microtubules.
The present paper is organized as follows.In Sections 2 we briefly introduce the mathematical model.More details of the biological background as well as its discussion may be found in [1].Section 4 is devoted to the mathematical analysis of the model.Finally we discuss about further improvements in Section 5.

Intracellular dynamics
Here we present a very brief description of the underlying biological process.If interested in more details, please refer to our work Szymanska at al. [1] or, to a textbook of molecular biology by Alberts at al. [12].
In simplest terms, an eukaryotic cell consists of a nucleus and cytoplasm, and the nucleus is located inside the cytoplasm, see Fig. 1.It means that considered domain, i.e. cell volume denoted by Ω, is composed of two separated parts -nucleus denoted by N and cytoplasm denoted by C .The entire cell is limited by a so-called cellular membrane here denoted by Γ c .The nucleus is separated from the cytoplasm by a membrane, which is called the nuclear envelope that is denoted by Γ n .Both on the cell membrane and the nuclear envelope we assume boundary conditions.In the cytoplasm there are so called microtubules, which are thin polymers arranged concentrically from the centre to the edges of the cell.They are used to give shape to the cell, and to facilitate intracellular transport of molecules and organelles.In the model the microtubules density is described by ξ function.The nucleus contains DNA, which is used for the production of mRNA molecules.The newly produced mRNA molecules escape from the nucleus to the cytoplasm.After passing the nuclear envelope, with help of the microtubules, they are moving into the cytoplasm, where they are released.Thus, in the present paper the mRNA in the cytoplasm is divided into associated with microtubules and free denoted by m i and m i , i = 1, . . ., N, respectively.In the cytoplasm there are many molecular machineries called ribosomes where, on the basis of mRNA, proteins are produced.Newly synthesised proteins diffuse in the cytoplasm where they perform its appropriate functions.Some of the proteins are needed in the nucleus in order to regulate the process of mRNA production.In such a case they bind to the microtubules and are transported along them to the nucleus.As in the case of mRNA, the protein can be divided into the associated with microtubules denoted by p i and free, denoted by p i , i = 1, . . ., N. In the model, the protein produced in ribosomes in the cytoplasm acts also as its own inhibitor.That means protein molecules are transported into the nucleus where they act as a suppressor for its own mRNA production.Concentration of mRNA in the nucleus is denoted by m i and concentration of protein there is denoted by p i , i = 1, . . ., N. This mechanism is called negative feedback and allows to maintain the limited protein concentration.

The system of equations
The cell volume, defined by Ω is a nonempty convex bounded and smooth domain in R d with d ≥ 1.In what follows we do not impose any restriction on the size of the space dimension d.The boundary, denoted by Γ c , corresponds to the cell membrane and c is the unit normal vector outward the cell.We defined the nucleus of the cell by N , a nonempty bounded and smooth sub-domain in Ω.Its boundary Γ n corresponds to the nuclear envelope, n is the unit normal vector outward the nucleus and we assume that Γ n ∩ Γ c = ∅.The cytoplasm C is defined as the complement of the nucleus in the cell, i.e.C = Ω \ N .Here and in what follows we adhere to the convention that given set U we denote by U T = [0, T ] × U, and T is arbitrary but fixed.In addition, the cytoplasm of eukaryote cells presents numerous microtubules which are thin filaments going from the nucleus envelope to the cell membrane.To model the microtubules, we introduce the unit vector η : C → R d oriented from the nucleus envelope to the cell membrane such that the microtubules are tangent to η.From now on, we denote by t > 0 the time and by x ∈ Ω the position in the cell.
In the present paper we consider the system of nonlinear PDEs with specific boundary conditions that can be referred to specific signaling pathway composed by N couples of mRNA and protein.The main novelty of the Γ n -nuclear envelope; ξ -concentration of microtubules; m i -mRNA molecules in the nucleus; p i -protein molecules in the nucleus; m i -mRNA molecules in the cytoplasm that are linked to the microtubules; m i -mRNA molecules in the cytoplasm that are free, i.e. not linked to the microtubules; p i -protein molecules in the cytoplasm that are free, i.e. not linked to the microtubules, i = 1, . . ., N; model presented in [1] comes from the distinction of two ways of motion in the cytoplasm.We propose to model the distribution of microtubules through the function ξ on C and such that The density of microtubules is supposed to be given at any time thus it is not an unknown of our model.The movement of molecules in the cytoplasm depends on whether the molecules are bound to the microtubules or not.Therefore, in the cytoplasm, both species, protein and mRNA are divided into those associated with microtubules, namely linked molecules, and those not associated, namely free molecules.For a given t and x we denote by m i (t, x) and p i (t, x) the concentration of mRNA linked to the microtubules and the concentration of protein linked to the microtubules, respectively.Analogously, by m i (t, x) and p i (t, x) we denote the concentration of mRNA and the concentration of protein, respectively, that are free to the microtubules.Thus, we consider the following equations in the cytoplasm with the nonlinear parameters that are local functions of the unknown concentrations and the density of microtubules.More precisely, the diffusion parameters k i (ξ; m c ; p c ) and κ i (ξ; m c ; p c ) are square matrices, the speed parameters v i (ξ) and ν i (ξ), the parameters i (ξ; m c ; p c ) and λ i (ξ; m c ; p c ), the exchange parameters b i (ξ; m c ; p c ) and β i (ξ; m c ; p c ), the transcription rate a i (ξ; m c ; p c ) and the translation rate α i (ξ; m c ; p c ) are positive scalar function of the vectorial unknowns In the nucleus, i.e. if x ∈ N , the concentrations of the i th couple are denoted by m i (t, x) for mNRA and p i (t, x) for protein.We assume that the motion in the nucleus is govern by the diffusion and the evolution of the concentrations are given by the following reaction-diffusion equations with the nonlinear parameters that are local functions of the unknown concentrations.More precisely, the diffusion operators k 0i (m n ; p n ) and κ 0i (m n ; p n ) are square matrices, the decay parameters 0i (m n ; p n ) and λ 0i (m n ; p n ) and the transcription rate a 0i (m n ; p n ) are positive scalar function of the vectorial unknowns m n = (m i ) 1≤i≤N and p n = (p i ) 1≤i≤N .At the nuclear envelope, i.e. if x ∈ Γ n , we suppose that the exchange is only one-way.More precisely the mRNA concentration m i leaving the nucleus is split into the subpopulations with concentrations m i and m i accordingly to the density of microtubules ξ at the nuclear membrane.The equations read as The protein concentration in the nucleus p i results from the contribution of two ways of motion, i.e.
with the speed passing through the nuclear nuclear envelope of mRNA u i (ξ) and of protein µ i (ξ).
Finally we consider no exchange of particles between the cell and its environment.It leads to the following boundary condition at the cell membrane Note that the boundary conditions (4) imply that the mRNA particles reaching the cell membrane via the advection are released into the cytoplasm and change its way of motion.We refer to [1] for details of the biological interpretation of Eqs. ( 1)-( 2) and the boundary conditions ( 3)-(4).

Analysis
In the following section, we give the main result of the paper.We prove the global well-posedness of the nonlinear and multi-dimensional model.The model is a system of non-linear coupled hyperbolic and parabolic equations and its well-posedness, especially globally in time, is not obvious.The proof is based on an adapted iterative approximation process which permits to use the classical results for the linear hyperbolic equations as well as those for the linear parabolic equations.The convergence of the approximations leading to the global well-posedness is obtained using classical contraction arguments.
Then, for all initial data such that there exists an unique global weak solution of (1)-( 2) with the boundary condition (3) and (4) such that Assumptions of Theorem 1 seem complex, in particular for the transcription rate a 0i , the translation rate α i and the decay parameters 0i , λ 0i , i and λ i .However, several models of signalling pathways with negative feedback proposed in the literature satisfies those conditions.For instance, this is the case of Hes1 gene regulatory network (Hes1 is a mammalian protein that suppress transcription) and p53/mdm2 regulatory network (p53 regulates the cell cycle and is one of the mayor tumour suppressors) [9].
In order to prove Theorem 1, we start with the following results: Let Ω ⊂ R d be a smooth bounded domain and be n the outward normal to its boundary ∂Ω.Assume that (L1h1) There exists ε > 0 such that for any (t, x) ∈ Ω T , the eigenvalues of the square matrix K (t, x) are larger than ε; (L1h2) For any (t, x) ∈ Ω T , we have M (t, x) ≥ 0; Then, for any initial data φ 0 ∈ L 2 (Ω), there exists an unique weak solution φ in L ∞ 0, T ; L 2 (Ω) ∩ L 2 0, T ; H 1 (Ω) of the following parabolic equation with the Robin boundary condition Proof.Let (ψ n ) n≥0 be a dense set in H 1 (Ω) with ψ 0 = φ 0 and H 1 n = Vect (ψ i ) 0≤i≤n .We consider the following approximation of the variational formulation of (5): For any initial condition in φ 0 ∈ H 1 (Ω), find φ n ∈ C 1 0, T ; H 1 n such that for any ψ ∈ H 1 n we have The approximated problem is a system of linear ordinary differential equations in the finite dimensional space H 1 n .We conclude that it is well-posed for any time T > 0. By setting ψ = φ n , using (L1h2), the Cauchy-Schwarz's inequality, the trace theorem and the binomial expansion, we obtain the following a priori estimate where C is the constant of the trace theorem, i.e.
Then there exists a subsequence weakly converging in L 2 0, T ; H 1 (Ω) and weakly-in L ∞ 0, T ; L 2 (Ω) to the weak solution φ of the initial problem (5) satisfying the a priori estimate (6).Since the problem is linear, the difference between two solutions satisfies the same problem ( 5) with a vanishing source term F = 0 and incoming flux G = 0. Using the a priori estimate (6), we conclude the uniqueness of the solution.
Lemma 2. Let Ω ⊂ R d be a smooth bounded domain and n be the outward normal to its boundary ∂Ω and Assume that (L2h1) For any time t ∈ [0, T ], we have V (t, .)∈ C 1 (Ω); Then, for any initial data φ 0 ∈ L 2 (Ω), there exists an unique weak solution φ in L ∞ 0, T ; L 2 (Ω) ∩L 2 (∂Ω T ; |V • n| dσdt) of the following hyperbolic equation Proof.Since C 1 (Ω T ) is dense in L 2 (Ω T ), we can define the sequences F n ∈ C 0 Ω T and G n ∈ C 1 ∂Ω − T ), bounded in L 2 (Ω T ) and converging to F in L 2 (Ω T ) and G in L 2 ∂Ω − T , respectively.Then for any initial condition To prove existence, it is enough to see that the following function is a solution with the flow χ (s, t, x) defined as the solution of The flow χ (s, t, x) is well-defined accordingly to the Picard-Lindelöf theorem and (L2h1).The time 0 ≤ τ (t, x) ≤ t is defined such that χ (τ (t, x) , t, x) ∈ ∂Ω − or τ (t, x) = 0. Then we define a sequence of initial data φ 0 n ∈ C 1 Ω , bounded in L 2 (Ω) and converging to φ 0 in L 2 (Ω).Multiplying Eq. ( 8) by the solution φ n we get the a priori estimate Neglecting the second term of the LHS and using the Grönwall lemma, (L2h1) and (L2h2), we conclude that the sequence φ n is uniformly bounded in L ∞ 0, T ; L 2 (Ω) .Then, there exists a subsequence weakly converging in L ∞ 0, T ; L 2 (Ω) to the solution of (7).Since the problem is linear, the difference between two solutions satisfies the same problem (7) with the vanishing source term F = 0 and incoming flux G = 0. We conclude therefore uniqueness of the solution φ of (7).Finally using the a priori estimate ( 9), we also get the estimation of the solution at the outgoing boundary ∂Ω + T .
Proof of Theorem 1.In the first step of the proof we use Lemma 1 and Lemma 2 to show the existence and uniqueness of the solution of the linearized model of ( 1)-( 2)-( 3)-( 4), more precisely assuming that the parameters are not function of the unknown functions.In the second step, we use contraction argument to show the convergence of the sequence of successive approximations and conclude the global existence of solution of the nonlinear system ( 1)-( 2)-( 3)-( 4).
Well-posedness of the linearized system: In this step we introduce successive approximations (m j i ), ( m j i ), (m j i ), (p j i ), ( p j i ) and (p j i ) of the solution of the linearized system, obtained by taking the source and decay terms with the previous approximations, and the boundary condition with the new approximation.More precisely, we take Assume that the previous approximations are in the corresponding Banach space Note that accordingly to the assumptions of Theorem, the initial data, considered as constant functions of time, are in the corresponding Banach space.We define the next approximation m j+1 i by We are going to show that the RHS is in L 2 (N T ).We note that in particular the estimate of term j 0i m j i = 0i m j n , p j n m j i is not trivial.In order to this end we use the Lipschitz continuity of the function 0i (m n , p n ) m i (assumed by (T1h6)).
However, the assumption that the function 0i (m n , p n ) m i vanishes when the unknowns vanish, i.e. 0 = 0i (0, 0) 0 = 0i 0, given by (T1h5), is needed.Then we obtain Similarly, the term a j 0i = a 0i m j n , p j n is estimated adding the term when the unknowns vanish, i.e. a 0i (0, 0) = a 0i .We have This yields the following estimate of the RHS This is clearly integrable since the previous approximations m j i and p j i are in L 2 (N T ).It leads that the auxiliary problem (10) satisfies Assumption (L1h3).We can easily check that Eq. ( 10) satisfies the other assumptions of Lemma 1.We conclude that there exists the unique global solution Then, we consider the approximations given by Similarly as for the auxiliary problem (10), we show that the RHS is in L 2 (C T ).In addition since m j+1 i is in L 2 0, T ; H 1 (N ) and using the Trace theorem we have T since u i is bounded (T1h3) and we conclude applying Lemma 2 that there exists the unique global solution Moreover, we consider Using the properties of m j+1 i and m j+1 i , by Lemma 1 we conclude that there exists the unique global solution Similarly we define the approximations by Using Lemma 1, we conclude that there exists the unique global solution Then we consider and by Lemma 2 we conclude that there exists the unique global solution Finally we consider and by Lemma 1 we conclude that there exists the unique global solution Global Well-posedness: This step is the proof of convergence of the sequences (m j i ), ( m j i ), (m j i ), (p j i ), ( p j i ) and (p j i ) when j tends to infinity.In the following, we use the sign function defined in a standard way In addition we denote We estimate the above differences by multiplying each equation by the sign of the difference and integrating over the space domain.More precisely, we multiply Eq. ( 10) by sgn δm j+1 i and by Lipschitz-continuity of the RHS (T1h6), we obtain Next we multiply (11) by sgn δm j+1 i and we have Multiplying ( 12) by sgn δ m j+1 i yields Summing up these estimates and using the boundaries balance we obtain Similarly we estimate the time evolution of the differences δp j+1 i , δp j+1 i , δ p j+1 i .Then summing with respect to 1 Finally by induction we obtain Since the RHS tends to zero when j goes to infinity, the sequences (m j i ), ( m j i ), (m j i ), (p j i ), ( p j i ) and (p j i ) converge.Then the global well-posedness of the nonlinear system (1)-( 2)-( 3)-( 4) follows.
We argue now the biological relevance of the model by the following positivity and conservation results: Proposition 1. [Global mass laws] The global concentration of each species satisfy the following equation with the total concentration of mRNA defined by and the total concentration of proteins defined by Proof.The result follows by direct estimates by the boundaries (3) and ( 4) balance.
The conservation law presented in Proposition 1 shows that the only processes that change of the global concentration of mRNA or protein are the transcription, the translation and the decay of molecules.Assuming the space homogeneous case, the model can be reduce to a ODEs system, see [13,14,15].However, the space distributions of the molecules are not homogeneous since the productions are located in a part of the cell.We performed the numerical simulations with parameters values relevant from biological point of view, that illustrate the heterogeneity of the space distribution, see [1].
Proof.We use the positive and negative part function 2φ ± = |φ| ± φ.First we multiply (2a) by the negative part of the solution m − i and we integrate over the space domain.It follows By the construction the negative part function is nonnegative.Since the transcription rate is nonnegative, it leads that the right hand side is nonpositive Since the initial data are nonnegative, we conclude that m i is nonnegative almost everywhere.In addition, the negative part of the outgoing flux vanishes Then we multiply (1b) by the negative part of the solution m i − and we integrate over the space domain.We have Using the Cauchy-Schwarz inequality we obtain which vanishes accordingly to (27).The last term of the RHS of (28) is clearly nonpositive, thus Using the Grönwall lemma, and since the initial data are nonnegative, we conclude that m i are nonnegative almost everywhere, i.e. m − i L 2 (C ) = 0 .
In addition, using the estimate (28), we conclude that the negative part of the outgoing flux and the negative part of the exchange term vanish, i.e.
Finally, we multiply (1a) by m − i and integrate over the space domain.We obtain We show that the positive part of the RHS vanishes using Cauchy-Schwarz inequality and the previous estimates (27),( 29) and (30).It follows L 2 (C ) ≤ 0 , and since the initial data are nonnegative, we conclude that m i is nonnegative almost everywhere, i.e. m − i L 2 (C ) = 0. Similarly we show that p i and p i are nonnegative almost everywhere.
The non-negativity of the concentrations is clearly essential in the term of biological interpretation.The fact that the solution of the model can not be negative allowed to simulate and analyse extreme situation for example when some parameters tend to vanish or become very high.

Conclusion
In this paper we have presented the preliminary mathematical analysis of the non-linear, multidimensional and multi-species model proposed by Szymańska at al. in 2014 [1], including active transport along the microtubules.The main properties such as well-posedness, positivity and global evolution are highlighted.However there are still many open mathematical problems including the analysis with less regular parameters, in particular transcription rate or speed along the microtubules.This questions seem relevant in the biological point of view and required a trickier analysis.In addition, the existence of periodic solutions that has been observed numerically [1] are still an open problem and are an important issus on the understanding of cell homeostasis and many disease like cancer.

Figure 1 :
Figure 1: Scheme of the space decomposition of the cell.Designation: Ω -eukaryotic cell; C -cytoplasm; N -nucleus; Γ c -cell membrane;Γ n -nuclear envelope; ξ -concentration of microtubules; m i -mRNA molecules in the nucleus; p i -protein molecules in the nucleus; m i -mRNA molecules in the cytoplasm that are linked to the microtubules; m i -mRNA molecules in the cytoplasm that are free, i.e. not linked to the microtubules; p i -protein molecules in the cytoplasm that are free, i.e. not linked to the microtubules, i = 1, . . ., N;