Reaction–diffusion systems for spatio-temporal intracellular protein networks: A beginner's guide with two examples

Spatio-temporal dynamics of a variety of proteins is, among other things, regulated by post-translational modifications of these proteins. Such modifications can thus influence stability and biochemical activities of the proteins, activity and stability of their upstream targets within specific signalling pathways. Commonly used mathematical tools for such protein–protein (and/or protein-mRNA) interactions in single cells, namely, Michaelis–Menten and Hill kinetics, yielding a system of ordinary differential equations, are extended here into (non-linear) partial differential equations by taking into account a more realistic spatial representation of the environment where these reactions occur. In the modelling framework under consideration, all interactions occur in a cell divided into two compartments, the nucleus and the cytoplasm, connected by the semipermeable nuclear membrane and bounded by the impermeable cell membrane. Passive transport mechanism, modelled by the so-called Kedem–Katchalsky boundary conditions, is used here to represent migration of species throughout the nuclear membrane. Nonlinear systems of partial differential equations are solved by the semi-implicit Rothe method. Examples of two spatial oscillators are shown. Namely, these are the circadian rhythm for concentration of the FRQ protein in Neurospora crassa and oscillatory dynamics observed in the activation and regulation of the p53 protein following DNA damage in mammalian cells.

may be also different from cell to cell; they can be influenced by extracellular factors such as light, heat, abundance of stress agents and duration of their action, growth conditions, available resources and many others. On a molecular basis, roles of proteins in specific networks are influenced by protein-protein interactions either through posttranslational modifications (by attaching a phosphate or acetyl group, ubiquitin to a protein, etc.), or by various compound formation. An interesting example is the signalling network of the protein p53 that may elicit life and death decisions in a cell. Activity of p53 in such situations (meaning not only transcriptional activity towards pro-arrest, proapoptotic and pro-survival target genes) is closely controlled by posttranslational modifications (such as phosphorylation or ubiquitination) and interactions with other proteins [1][2][3].
From a mathematical point of view, post-translational modifications can be advantageously modelled as enzyme reactions where a modifying protein (providing a phosphate group or ubiquitin) is considered to act as enzyme, i.e. it binds a substrate to form a compound, modifies the substrate in the compound and is released from the compound, becoming available to modify other target substrates. Mathematically, enzyme reactions can be described by using the law of mass action and the quasi-steady-state approximation yielding thus a non-linear term for a gain of the modified protein (the product of the reaction) and a loss of the former state of the protein (the substrate entering the reaction) [4,5]. This approach turns enzyme reaction kinetics into ordinary differential equations (ODEs) broadly used by mathematicians.
One of the simplest 'protein-mRNA' models is a non-linear biochemical oscillator proposed by Goodwin in 1965 [6] that simulates expression of a single gene controlled by its protein product P, thus, closing a negative feedback loop. In particular, equations appearing in the original Goodwin model are, respectively, reporting here only production (transcription in the first equation and translation in the second) and degradation terms (δ 1 and δ 2 ). Based on this model, a compartmental ODE model for the circadian clock rhythm of the FRQ protein expressed during a day in Neurospora crassa was proposed by Leloup and Goldbeter [7]. The model includes expression of the frq gene in the nucleus, a process that is influenced by available light, and translation of FRQ mRNA into FRQ in the cytoplasm that, afterwards, enters the nucleus and negatively regulates the transcription of its gene, as it is observed in vivo in [8].
More complicated 'protein-protein' and 'protein-mRNA' ODE models have been designed to represent intracellular signalling of the protein p53 [9][10][11][12][13][14][15][16][17][18]. One of them is also proposed in [19] and it simulates initial activation and regulation of p53 in response to DNA damage. Activation is performed through the interactions with the ATM protein, that is identified as DNA double strand break (DSB) sensor [20], and regulation of p53 by the proteins Wip1 and Mdm2, on whose genes p53 act as a transcription factor [21,22]. These four proteins have been identified as necessary to create a minimal p53 network yielding sustained oscillations in p53 concentration in vivo [9,23].
The spatial organisation of the cell, however, suggests to consider migration of the involved species in the cytosol and between the compartments, and thus to model a protein signalling more realistically by including diffusivities of the species and by setting particular translocation conditions to model exchanges of the species between the compartments. Thus, one has to deal with coupled systems of nonlinear evolution partial differential equations (PDEs) for the concentrations of the species in the nucleus and in the cytoplasm with transmission boundary conditions (BCs) imposed on the inner nuclear envelope and on the outer cellular membrane (Robin-like and zero-flux BCs in our modelling setting). Among other mathematical methods that are available in the literature, the semi-implicit Rothe method [24] can be used when dealing with non-linear systems of evolution equations. Without going into details, let us mention that after a suitable linearisation of non-linear equations, one can discretise time derivatives and solve a finite number of elliptic equations for which a non-negative and unique solution is guaranteed by the Lax-Milgram theorem [24,25]. A Rothe function formed from the solutions of the elliptic problems then can serve as a good approximation of a solution of the original reactiondiffusion problem (that is also non-negative and unique) in an entire time interval for the problem under consideration.
The aim of this article is to acquaint readers with reaction-diffusion equations for protein spatio-temporal signalling rising from proteinprotein (and/or protein-mRNA) interactions in the two compartments of a cell, nucleus and cytoplasm, and from migration of the species in and between these compartments. Since it is impossible in general to derive analytical solutions for such coupled systems, numerically, passing to a weaker notion of solution, it can be proved to exist even in Lipschitz domains (as cells are assumed to be in our approach). With the Rothe method in hand, examples of such spatio-temporal oscillatory models, particularly, the Leloup-Goldbeter FRQ model and ATM-p53-Wip1-Mdm2 model in individual cells are shown, numerically solved and illustrated.
Let us mention here that works on reaction-diffusion models applied to intracellular biology had been rather scarce until recently. Among previous studies [11,[26][27][28][29], it is worth pointing out the Ph D thesis of A. Serafini [30] that gives a different (theoretical and computational) approach to the reaction-diffusion equations applicable in modelling protein networks.
The organisation of the paper is as follows. Kinetics of enzyme reactions is briefly summarised in the following Section 2. Michaelis-Menten and Hill kinetics are presented as tools for mathematical modelling. An introduction to reaction-diffusion equations is given in Section 3, followed by two sections where the Leloup-Goldbeter model (Section 3.1) and the p53 model (Section 3.2) are demonstratively shown. Note that detailed analysis of the models is not the aim of the article. Whereas, to our best knowledge, the reaction-diffusion Leloup-Goldbeter model is newly presented here, so that we believe that no analysis has been done so far, an analysis of the p53 model has been given in [31]. Moreover, we describe in Section 3.3 the semi-implicit Rothe method that can be useful in numerical simulations. Finally, a short discussion closes the article.

A brief tutorial on the kinetics of enzyme reactions
Let us briefly summarise basic mathematical ideas used to model an enzyme reaction where a substrate S reacts with an enzyme En forming a complex SEn with a rate k 1 . In this complex, the enzyme En converts S into a product protein P with a kinetic rate k 2 . The enzyme En is released and available for further reactions. The first reaction between S and En is reversible, meaning that the complex SEn can eventually fall apart with a reverse rate of reaction k −1 , whilst the second reaction is assumed to proceed in one direction only [5,32]. The law of mass action for the concentrations of species denoted subsequently by s = [S],e = [En],c = [SEn] and p = [P] gives four equations, one for each species. These are and we can assume initial conditions s(0) = s 0 , e(0) = e 0 , c(0) = 0 and p(0) = 0. Note that the last equation can be easily solved whenever c is determined. In this case p(t) = k 2 ∫ 0 t c(τ)dτ and so the last equation can be omitted from further considerations. Note also that de dt þ dc dt ¼ 0 and so e + c = e 0 where e 0 is the total amount of available enzyme. By substituting e = e 0 − c into the equations for s and c, we can write The quasi-steady-state approximation (QSSA) can be finally applied to eliminate the equation for the complex. In particular, one can approximate the rate of change of c by zero, i.e. dc dt ≈0. Although equality is not always true during the reaction, see for example [4], by taking the right hand side of Eq. (2) to be zero, one can write c in terms of s and such term substitution into Eq. (1) finally yields where k 2 is the kinetic rate of the reaction and K 2 is affinity constant. A basic assumption for the QSSA is that the concentration of enzyme is smaller than the concentration of substrate (e.g., for ε = e 0 /s 0 small, typically, between 10 − 2 and 10 − 7 , [5]), the assumption that is not usually true in our models since enzymes under consideration may be present in large quantities. The QSSA can however be found valid also in the case when the concentration of the intermediate complex SEn is immediately consumed, i.e. it is always small compared with the total concentration of the substrate [4,5]. For a rigorous mathematical treatment of the QSSA from the point of view of perturbation theory, see [4].
Eq. (3) is thought to represent the total loss of the substrate S in the reaction, and thus, by using conservation of mass, the same term with opposite sign can be written as the total gain of the product P. In practical simulations, e 0 in Eq. (3) is often replaced by the actual concentration e. Example 2.1. In a model discussed in Section 3.2, we will consider the protein p53, its activation by the kinase ATM p (phosphorylation of p53 by ATM giving p53 p ) and regulated by the phosphatase Wip1 and the E3 ligase Mdm2 (dephosphorylation of p53 p by Wip1 and ubiquitination of p53 by Mdm2, respectively). All these modifications can be represented as enzyme reactions, i.e.
Complex → k ub Mdm2 þ p53 ub ubiquitination of p53 by Mdm2 ð Þ : Following the approach described above, one finally arrives at the equations for p53 and its phosphorylated version p53 p , in particular, with specific kinetic maximum velocity rates k x and affinity constants K x . The situation is slightly complicated whenever monomerisation and polymerisations of a protein are modelled as enzyme reactions, i.e.
where a polymer S consisting of N molecules reacts with an enzyme En yielding N monomers P in the first reaction, and N monomers S react with an enzyme En to give a polymer P in the second equation. In the same p53 model in Section 3.2 we will consider ATM dimer/monomerisation, since ATM in inactive state preferentially forms dimers and these dimeric molecules dissociate into active ATM p molecules once they sense occurrence of DNA damage [20], see also [19] for the derivation of equations for ATM dimer/mono-merisation. Kinetics of enzyme reactions using the law of mass action and the QSSA is sometimes collectively called Michaelis-Menten kinetics [4,5].
In the models presented below, we will also consider the expression of genes, particularly transcription of frq, mdm2 and wip1 into messenger RNA and its translation into proteins. Whilst translation is usually modelled as linear gain of the concentration of translated proteins, transcription can be modelled as enzyme reaction; however, transcription factors often express cooperative binding to the active sites of DNA (considered to be the enzyme in the reaction), i.e. binding of one substrate molecule affects binding of subsequent substrate molecules, as it is in the case of p53 that bind DNA as tetramers and thus four molecules of p53 cooperate whilst performing the transcription job [33,34]. In such kinetics, sometimes called non-Michaelis-Menten kinetics or Hill kinetics, the Hill equation can be used to describe the degree of cooperativity quantitatively [5,35]. In particular, transcription of mRNA can be described by a Hill function with Hill coefficient n, rate of transcription k t , affinity constant K t and a transcription factor TF. According to [34], see also [5] chap. 1, the Hill coefficient is equal to the number of binding sites, which, for example, in the case of tetrameric p53 binding DNA, is n = 4. Finally, natural degradation terms will be modelled here by linear loss terms or by Hill functions of coefficient 1. In the p53 model there is also a decay term in the nuclear and cytoplasmic p53 equation coming from ubiquitination of p53 by Mdm2, thus forming Michaelis-Menten kinetics discussed above. This is because p53 is degraded by the ubiquitin-dependent degradation machinery in many conditions [22]; for the sake of simplicity, we will not introduce any special equation for ubiquitinated p53 since it would consist only of degradation terms for this species.

Reaction-diffusion models for protein cellular signalling
Having in mind to represent the signalling of any protein in individual cells whose activity is determined by its post-translational modifications, one can write as many reactions as necessary to sufficiently describe a desired protein dynamics. Michaelis-Menten and Hill kinetics turn protein-protein interactions and transcription of genes into proteins into a system of ODEs. Since the spatial representation of signalling pathways can reveal diffusion patterns that are hidden in ODEs, one might consider a system of PDEs instead of ODEs. In PDEs, the protein dynamics in cells is particularly endowed with diffusivity rates and with some transmission conditions on the nuclear and cell (possibly other) membranes to represent, respectively, migration of species in and their exchange between the cellular compartments and between cells themselves. In the two sections that follow we show two examples of oscillators given by reaction-diffusion PDEs [36,37]. The first one is based on the commonly known and used Leloup-Goldbeter ODE model for circadian rhythms of the protein FRQ in Neurospora [7]. 1 The second one is an ODE and PDE model for the p53 intracellular dynamics following DNA damage in mammalian cells. Whilst the p53 ODE model has been proposed in [19], the PDE model, based on the ODE model, has been recently studied in [31]. Another example of a signal transduction modelled by reaction-diffusion equations, more precisely, a model of Ran-driven transport process was studied in [30].
In cellular settings, a general coupled reaction diffusion system may be written as follows, are the nuclear and cytoplasmic concentrations of N chemicals, T N 0, N N 0, f:ℝ N → ℝ N and g:ℝ N → ℝ N collect protein-protein and protein-mRNA reaction terms, and production/degradation terms for all N species as they are discussed in Section 2, D is a diagonal N-by-N matrix with the diffusivities on the diagonal, and div is the divergence operator. The domains Ω 1 and Ω 2 represent the nucleus and the cytoplasm, respectively, assumed to be Lipschitz domains [24,40], Γ 1 is the nuclear membrane and Γ 2 is the cell membrane. The geometry of a cell under consideration is sketched in Fig. 1, that in our very simplified framework represents a cell as big as 10-100 μm or so, which is able to proliferate and which can have heterogeneous shape. However, we exclude cells of complicated structures and morphologies such as neurons or muscle cells.
As previously outlined, we will consider zero-flux boundary conditions on the outer cell membrane Γ 2 , i.e. ∂vi ∂n2 ¼ 0 for i = 0,1…, N − 1 where n 2 is the unit normal vector oriented outward from the cell, Fig. 1. Although species with molecular weight over 40 kDa usually use active transport mechanisms to be translocated from one compartment to another, for the sake of simplicity, only passive transport process driven by the difference in concentrations at both sides of the nuclear membrane, modelled by the so-called Kedem-Katchalsky BCs, is used here [11,31,41,42]. In particular, export/import of the nuclear species is modelled by and export/import of the cytoplasmic species by where n 1 is the unit normal vector oriented outward from the nucleus (therefore, the BCs in Eq. (6) come with a different sign), and p i are the permeabilities of the membrane for each species i = 0,1…, N − 1.
Recall that if a chemical is assumed to migrate between the compartments in one direction only, e.g. mRNA is supposed to move from the transcription sites in the nucleus to the cytoplasm and not the other way round, then this can easily be modelled by omitting either v i or u i in Eqs. (5) and (6), cf. (9)-(11) and (14) below. In addition to the boundary conditions, we will assume initial conditions for each species i = 0,1…, N − 1 to be non-negative functions 3.1. The Leloup-Goldbeter model for the N. crassa circadian clock 3.1.1. The frequency protein and circadian clock Circadian clock in Neurospora controls the pattern of asexual development [43] and a central component of sustained rhythms is the FRQ protein (a long form of 989 amino acids and a short form of 890 amino acids) encoded by the frq clock gene [8]. The oscillator consists of a self-control of frq transcription, i.e. a negative feedback loop in which frq is synthesised into the FRQ protein (both forms are required for robust rhythmic responses) which, in turn, acts to reduce intensity of frq transcription into mRNA [43].
Making the Neurospora clock start at midnight when FRQ mRNA and FRQ proteins attain low levels in concentration, we see that transcription of the gene increases and mRNA reaches its peak in the midmorning about 4 h before the peak of the total FRQ. FRQ enters the nucleus where it contributes to inhibit gene transcription through the interaction with other factors (WC-1 and WC-2) required in the Cell scheme: the nucleus Ω 1 , the cytoplasm Ω 2 , the nuclear membrane Γ 1 and the cell membrane Γ 2 ; n 1 and n 2 are the unit normal vectors oriented outward from Ω 1 and Ω 2 , respectively. Table 1 Parameters for the FRQ PDE model. Note that the diffusivity of the mRNA-protein complex is D 0 = 108 μm 2 h −1 as suggested in [48] and the diffusivity of FRQ is D 1 = 43,200 μm 2 h −1 that corresponds to the diffusivity of a protein with the weight 110 kDa [41,31], whilst FRQ has the molecular weight in the range 97-160 kDa depending on its form [44,49]. The permeabilities p 0 = 1 μm h −1 and p 1 = 2.7 μm h −1 are tuned by hand so that the system gives oscillations.

Parameter
Value

Table 2
Parameters for the p53 PDE model. Note that the diffusivity of the mRNA-protein complexes is D 2 = D 6 = 1.8 μm 2 min −1 as the diffusivity of an average mRNA-protein complex is estimated to be in range 1.2-2.4 μm 2 min −1 in [48]. The diffusion rate of p53 is chosen to be D 0 = D 3 = 1000 μm 2 min −1 which is in agreement with the estimated value for p53-GFP in [55]. Other diffusivities are chosen with respect to the weights of the species corresponding to the weight of p53. More details about the parameter choice are provided in [31].
frq gene transcription [43,44]. Partial phosphorylation of FRQ leads to its initial degradation as soon as it enters the nucleus. Degradation then increases with more extensive FRQ phosphorylation through the early night [43,45].

Spatio-temporal Leloup-Goldbeter PDE model for FRQ
Circadian clock represented by sustainably oscillating concentration of the FRQ protein is modelled by Leloup and Goldbeter [7] via the following ODE equations, The first term in the first equation in Eq. (7) models the transcription of the frq gene that is negatively controlled (in a steep way if the Hill coefficient r is high) by the expression level of the nuclear FRQ protein, and where the rate of gene transcription V s is a function of light available during the day and changes in time, V s = V s (t,light). The second term represents mRNA degradation with rate V m . The first term in the second equation stands for gene translation into the cytoplasmic protein, and it is followed by the degradation term and the terms modelling exchange of the protein between the compartments. The equation for the nuclear concentration of the protein contains only terms for its exchange between the compartments.
The model of Leloup and Goldbeter [7] is based on the Goodwin model [6] and Hill functions are used to simulate production/ degradation events occurring in the cell. Starting with the model (7), we can formulate reaction-diffusion equations for the nuclear and cytoplasmic FRQ protein and its nuclear and cytoplasmic mRNA. To our best knowledge, this is the first attempt to model the classical Leloup-Goldbeter circadian oscillator [7] via partial differential equations. For simplicity, let us simplify notations and write for the nuclear and cytoplasmic concentrations of FRQ and mRNA of FRQ. The PDE equations may become where, compared with Eq. (7), degradation of FRQ and its mRNA is considered to occur in both compartments with the (possibly different) rates of degradation V m n , V m c , V d n , and V d c (and the corresponding affinity constants) in both compartments distinguished by the superscripts n and c for the nucleus and the cytoplasm, respectively. Note that, similarly as we modelled the translation process in [31], and as it was also modelled in the spatio-temporal Hes1 model in [26], translation of mRNA is assumed to occur at distance from the nucleus. More precisely, it is observed that the proteins that move from the translation sites in the cytoplasm back into the nucleus are likely translated outside of the endoplasmic reticulum (ER) [26,46]. Translation of FRQ mRNA into proteins is thus assumed to occur in a subdomain of the cytoplasm denoted by χ tra ; cf. Fig. 4(a) in Section 3.3.
Migration of the FRQ protein can be modelled in both directions from the nucleus to the cytoplasm, and the other way round, by the Kedem-Katchalsky BC The flux in Eq. (9) is determined by the difference between the nuclear and cytoplasmic concentrations at the nuclear membrane, to be compared with translocation in the ODE model in Eq. (7), 2nd and 3rd equation, that runs in both directions with generally different rates k 1 and k 2 . In addition, there is also a biological evidence [47] that FRQ is localised mainly in the nucleus and so we can also assume that it can only move from the translation sites in the cytoplasm (ribosomes) into the nucleus. In this case the Kedem-Katchalsky BC is The mRNA is assumed to move from the nucleus to the cytoplasm where it binds free ribosomes localised outside of the ER, thus the boundary conditions on the nuclear membrane are set to where n 1 is the unit normal vector pointing outward from the nucleus and p 0 and p 1 are the permeabilities for FRQ mRNA and its protein product, FRQ. Zero-flux boundary conditions on the cell membrane ∂v i ∂n 2 ¼ 0, i = 0,1 are still assumed.
Nuclear and cytoplasmic concentrations as solutions to the PDE model (8) with (10) and (11) are shown in Fig. 2 for V s set to 1.6 nM h − 1 in the continuous darkness, Fig. 2(a) and for V s defined as a piecewise constant function attaining values 2 nM h − 1 and 1.6 nM h − 1 every 12 h in the 12:12 light:dark simulations, Fig. 2(c). Periods of oscillations are 21.5 h and~24 h in constant darkness and 12:12 light:dark simulations, that can be compared with the periods of oscillations in the Leloup-Goldbeter ODE model and experimentally observed periods in Neurospora, [7] and citations therein. Fig. 2(b) and (d) shows stable limit cycles reached in the FRQ and its mRNA concentrations.

The protein p53 and its intracellular signalling in genome protection
The protein p53, the so-called guardian of the genome, plays important roles in genome protection since it is able to either trigger apoptosis or stimulate permanent cell cycle arrest; thus it protects the cell integrity from turning to malignancy, as well as it contributes to cell survival through initiation of some DNA repair processes; it thus contributes to both cell renewal and tissue repair. In contrast to the positive role of p53 in protection against cancer, p53 may be responsible for many unwanted effects in ageing, as well as in the debilitating toxic side effects of chemotherapeutic treatments [2,50]. The protein p53 thus acts as a hub in a broad range of various signalling pathways. Interested readers who might like to become acquainted with other reviews about p53 protein signalling are referred to [3,22,51,52].
The intracellular dynamics of p53, namely its activation and regulation in response to DNA damage caused by γ-radiation or some drugs in chemotherapy, involves the kinase ATM as a damage sensor and activator of p53, and the downstream p53 targets Wip1 and Mdm2 that regulate p53 through protein-protein interactions (phosphorylation and ubiquitination) and its subsequent degradation.
ATM in inactive state (preferentially) forms dimers that promptly dissociate into active monomers following occurrence of DNA DSBs [20] (modelled here as an interaction of dimeric ATM with an abstract signal E standing for the abundance of DNA damage [31]). Such mono-mers then phosphorylate p53 on serine 15 (Ser15) residue that is located very close to the Mdm2 binding domain [22]. Ubiquitindependent degradation controlled by the Mdm2 ligase is a principal way of p53 regulation; p53 molecules tagged by ubiquitin are exported from the nucleus to the cytoplasm and degraded [22]. Hence, phosphorylation of p53 by ATM impedes Mdm2 to bind p53 and so Mdm2 cannot ubiquitinate p53. Halted p53 ubiquitination results in p53 stabilisation in the nucleus where it forms tetrameric compounds that bind DNA and act transcriptionally towards many target genes, with the mdm2 and wip1 genes among them [33,34]. The wip1 gene protein product, the phosphatase Wip1, then dephosphorylates p53 on Ser15, thus it unmasks p53 from Mdm2. Wip1 also dephosphorylates ATM molecules that bind another dephosphorylated monomers to create inactive dimers [21,53,54]. All these four proteins have been observed to generate sustained oscillations in response to DNA DSBs in vivo, with the period of p53 oscillations in the range 4-7 h [9,12,23].

Spatio-temporal PDE model for p53
In a very simplified framework, intracellular activation and regulation of p53 thus consist of the two negative feedback loops, ATM-p53-Wip1 and p53-Mdm2. Representing protein-protein interaction and gene transcription by Michaelis-Menten and Hill kinetics discussed in Section 2 together with other assumptions on migration of the species throughout the nuclear membrane, one can formulate a reaction-diffusion system for the nuclear and cytoplasmic concentrations of the proteins/mRNAs [31], cf. Eqs. (12) and (13) below. Before stating equations for the PDE model, let us also show the ODE model for p53 that has been studied in [19]. Note that unlike the model studied in [19], the transcription of the mdm2 and wip1 genes is fully controlled by p53; thus we do not consider any constant rates for the Mdm2 and Wip1 mRNA production independent of p53, which were previously assumed in [19] and also in [11]. In both ODE and PDE models, natural degradation of the species and translation of mRNA into proteins are modelled as linear gain/loss and a constant term is used for p53 production, since, for sake of simplicity, we do not assume any promoters (transcription factors) for p53. Difference between the models is naturally in the diffusivities integrated into PDEs (which are supposed to be zero for each species in ODEs). Note that the Kedem-Katchalsky BCs, and also passive transport mechanism for exchange, cf. Eq. (14) below, can readily be rewritten into similar conditions in the ODE setting (with the permeabilities denoted similarly by p i for each species i).
Let us denote concentrations of the species by  Table 1. The plotted concentrations are scaled over the total volume of the computational domain (nucleus or cytoplasm).
The ODE system that, with minor changes, was previously studied in [19], is where the velocity of migration in the nucleus is balanced by the special volume ratio V r assuming that the cytoplasm has V r times bigger volume than the nucleus [10,19].
The dynamics of the species together with their spatial distribution in the nucleus can be described by the following seven equations, where, in contrast to the model in [31], the mdm2 and wip1 gene transcription is fully controlled by p53 here. Equations for the cytoplasmic  (12) and (13) the Kedem-Katchalsky BCs (14). In (c) and (d) stable limit cycles of the nuclear concentration of ATM p to the total concentration of p53, and the total concentration of p53 to Mdm2 are shown. Protein signalling is rescaled so that the period of p53 oscillations is 6 h. The parameters chosen in simulation are in Table 2, for their detailed description see [31,19] and also citations therein. The plotted concentrations are scaled over the total volume of the computational domain (nucleus or cytoplasm).
concentrations are rather simpler, since only degradation of the species and translation of mRNA, supposed to occur outside of the ER [26,46], are assumed to occur there. Free ribosomes needed for translation of Mdm2 and Wip1 mRNA are thus localised in a cytoplasmic subdomain denoted by χ tra . In addition, production of p53 (modelled by a constant term k S ) occurs in a subdomain identified by χ bp ; cf. Fig. 4(b) in Section 3.3. Cytoplasmic equations thus read The Kedem-Katchalsky BCs (6) set on the nuclear membrane Γ 1 are The PDE model (12) and (13) with the Kedem-Katchalsky BCs (14) has been previously proposed in [31]; we refer to this work for more information about p53 dynamics, for biological evidences that are behind the model assumptions, migration of the species throughout the cell. Modelling issues, parameter selection and results from simulations including analysis of the PDE model with respect to spatial perturbations and abundance of the DNA damage E can be found in [31] as well. Fig. 3 shows the nuclear and cytoplasmic concentrations of the proteins ATM, Mdm2, Wip1 and the total concentration of p53 (phosphorylated and unphosphorylated) with time rescaled so that the period of p53 pulses is 6 h.

The semi-implicit Rothe method for numerical applications
The semi-implicit Rothe method [56] can be used when one has to solve a coupled reaction-diffusion model (4) in cells as the ones previously discussed in Sections 3.1 and 3.2 as well as it can be used to prove theoretically existence of such solutions. Let T N 0 and I = [0,T] be a fixed time interval. Without going into details (a special article [40] is dedicated to more precise description of the method and rigorous proofs of statements claimed here), let us mention that we are looking for solutions u = [u 0 ,u 1 ,…,u N − 1 ] T and v = [v 0 ,v 1 ,…,v N − 1 ] T to the coupled reaction-diffusion system defined in two time-space cylinders I × Ω 1 and I × Ω 2 connected on a common boundary I × Γ 1 by Robin-like boundary conditions (Kedem-Katchalsky BCs) and satisfying zero-flux BCs on I × Γ 2 and some initial conditions at time t = 0. Recall that in d-dimensional case, d ≥ 1, Ω 1 and Ω 2 are open and bounded sets in ℝ d with Lipschitz boundaries Γ 1 = ∂Ω 1 ∩∂Ω 2 and Γ 2 = ∂Ω 2 − Γ 1 , see also Fig. 1 for a schematic representation of the domains. Solutions u i for i = 0,1,… N-1 are supposed to belong to a standard Sobolev-Bochner space, with V 1 = W 1,2 (Ω 1 ) and V 1 * dual to V 1 . Similarly, v i ∈W 2 1,2,2 (I;V 2 ,V 2 * ) with V 2 = W 1,2 (Ω 2 ) and V 2 * dual to V 2 for all i = 0,1,… N − 1. By du dt we understand the derivative of u in sense of distributions. We will refer to a solution of the coupled reaction-diffusion system as a solution pair and write w = {u,v} or In general, the reaction-diffusion Eq. (15) can be understood as an abstract initial-boundary value problem: find u(t) and v(t) solutions to Þ¼e g t ð Þ for a:a: t∈I where A i : of time t, and e f and e g independent of u and v, respectively; u 0 ∈L 2 (Ω 1 ) and v 0 ∈L 2 (Ω 2 ). Let us define τ N 0 a time step such that T/τ is an integer for simplicity. In the Rothe method [24,56], one has to discretise time t by backward differences and, if necessary, to approximate e f and e g at particular points t = kτ of a time grid (assumed to be equidistant for simplicity), k = 0,1,…,T/τ. This approximation can be done, for example, by convolution of e f and e g with suitable mollifiers [24,25]; however, reaction terms rising from protein-protein (protein-mRNA) interactions are smooth functions, thus well defined at each time point t.
To efficiently handle terms A and B, as they are non-linear in our reaction-diffusion systems, we can consider a certain linearisation A z; Á ð Þ : V 1 →V Ã 1 of A at a point z and B z; Á ð Þ : V 2 →V Ã 2 of B at a point z. Then, one can define u τ k and v τ k for k = 1,2,…,T/τ by where the conditions A u; u ð Þ¼A u ð Þ and B v; v ð Þ¼B v ð Þ are required. The initial conditions u τ 0 and v τ 0 are set to the original initial conditions u 0 and v 0 or to their suitable approximations. The Kedem-Katchalsky BCs on Γ 1 and the zero-flux BCs on Γ 2 for k = 1,2,…,T/ τ become can be used as approximations of solutions u and v to the original reaction-diffusion systems. More precisely, it can be shown that u(t) and vτ(t) weakly converge to u(t) and v(t) by limit passage for τ → 0 for a.a.t∈I. This non-rigorous procedure constructively proves existence of u(t) and v(t) in the time-space cylinders I × Ω 1 and I × Ω 2 , more details are left to [40]. Much richer knowledge, and not only on the Rothe method, can be found in the book of Roubíček [24].
Example 3.1. For convenience and because the Leloup-Goldbeter system is less demanding on space, we can show the equations of the model (8) with the Kedem-Katchalsky BCs (10) and (11) appearing in the Rothe method. In particular, the equations in Eq. (8) may be discretised and linearised into the following elliptic problems whilst the Kedem-Katchalsky BCs (10) and (11) become It is, however, worth mentioning that the computational time needed to obtain a solution to the Leloup-Goldbeter model by the semiimplicit Rothe method is 3.3-fold lower than the time required to solve the same problem with the Newton method (computations with T = 72 h and τ = 0.1 took 105 s in the semi-implicit Rothe method implementation and 344 s in the Newton method with the accuracy tolerance tol = 10 −3 needed to terminate iterations). Similarly, the Rothe method is 5.2-faster than the Newton method when computing solution to the p53 model (within t = 400 time units and τ = 0.2 took 37 min in the semi-implicit Rothe method implementation and 193 min in the Newton method with tol = 10 −3 ).
For the Leloup-Goldbeter and p53 reaction-diffusion systems presented in this article, the semi-implicit Rothe method and the Newton method have been implemented in the freely available solver FreeFem++ [57] and we ran 2D simulations on a disk-shaped cell with radius 10 μm, see Fig. 4 on a standard machine MacBook Pro, 2.4 GHz Intel Core icontain contain 7 processor and 8 GB (1600 MHz) memory.

Discussion and conclusion
In this short review we have raised the question of replacing commonly used ODE models for protein-protein interactions by reactiondiffusion models which are well studied [24,25,37] and which not only contain reaction terms coming from the protein interactions but also describe the spatial distribution of species involved in the reactions over the cellular compartments. We have proposed the simpler Leloup-Goldbeter circadian model for FRQ in Neurospora that contains two equations for the nuclear and cytoplasmic FRQ protein and its mRNA, and a more complicated model for p53 response to the occurrence of DSBs in a single cell. These models, and more generally reactiondiffusion PDE models for intracellular protein dynamics, are likely to be more realistic than ODE models.
As an example of such possibly more realistic models, let us consider the slightly more complicated PER model in Drosophila of Leloup and Goldbeter [38,39], that relies on the same principles as the FRQ model. Indeed, two additional phosphorylations are considered, and this first PER model (more sophisticated models of PER have been published by Leloup and Goldbeter) is proposed to be amenable to describe two Drosophila mutants for PER, with shorter or longer period of oscillations, by lower and higher values, respectively, of the maximum degradation rate v d of cytoplasmic a b  8) and (13). The nucleus, thus shown as an inner disk, has radius ffiffiffiffiffi ffi 10 p μm in both cases. The endoplasmic reticulum is an annulus with radii ffiffiffiffiffi ffi 10 p μm and 5 μm. The rest of the cytoplasm is an annulus with radii 5 μm and 10 μm. The basal production of p53 in the p53 model is assumed to occur in an annulus with radii 5 μm and 6 μm (inner peripheral, i.e., not perinucleic, zone, green) in (b), denoted by χ bp in 13. diphosphorylated PER. But is this sole change of constants the best reason, in a more physiological model, for the change in periods? Could it be related to transcription of genes, to species diffusion, to translation into protein products, to cytosolic/nucleic degradation, to nucleocytoplasmic transport, to the extension of the "dead zone" around the nucleus where translation occurs (Fig. 4a)? The same questions can be posed when one considers p53 and drugs amenable to electively modify specific parts of its intracellular dynamics. This can be rendered optimally in PDE models, that are much more physiological than (rougher) ODE models, since they are naturally able to take intracellular spatial features into account, including possible space heterogeneities in the intracellular medium (not considered in the examples, apart from the above mentioned dead zone of non-translation in the cytoplasm).
The oscillatory dynamics of p53, "the guardian of the genome", has long been evidenced by biological recordings in experimental conditions representing cell stress due to, e.g. radiotoxic insult [12,23]. In this sense, the biological question of identifying intracellular spatiotemporal dynamics is certainly crucial. Furthermore, since p53 disruptions are found in more than 50% of solid tumours, with various modifications of its stress-induced dynamics [2,3,50], it is also a crucial point in cancer therapy modelling to understand how p53 dynamics may be affected in disease, and to this aim, to have an accurate representation of how this dynamics works in physiological conditions, which is the purpose of designing p53 dynamic models.
We have shown how easy it is in principle to start from an ODE compartmental model, add diffusion termswhich in their simplest version are mere laplaciansand slightly modify the representation of exchanges between the compartments to adapt them to the new setting, to obtain a ready-made reaction-diffusion PDE model. However, dependence on spatial patterns, such as oscillations in concentration, on the diffusivities and the (nuclear membrane) permeabilities has to be better studied, with precise identification of the underlying biological parameters, to understand their roles properly. To obtain more realistic PDEs one can consider effects of different viscosities of buffers on different diffusivities of molecules or adjunction of advection terms when knowledge of aided transport makes it relevant (e.g., transport through microtubules). To support these ideas, different mobility of molecules in the nucleus and the cytoplasm caused by limited diffusions has led to biphasic caspase activation kinetics, see [58] and citations therein. In some circumstances it might be convenient to include specific spatial structures, such as position and density of ribosomes and thus to model gene expression more specifically [59,60].
In addition, intracellular signalling in our rather simplified modelling setting is restricted to cells of physiology and morphology where a molecular network in question is, at least, partially understood. What does the p53 signal transduction in response to DNA damage in a nerve cell or in a polynucleic muscle cell look like? What is the role of diffusion in a cell where a signal is spread over long distances and in very small cells. These and other questions should be further addressed and spatial PDE models might fruitfully be used for this purpose.
Note also that the oscillatory patterns are self-organised not only due to events occurring in the nucleus or in the cytoplasm but that they are also tightly connected to boundary conditions on the membranes. Thus physiological delays maintained due to the semipermeability of the nuclear membrane are more typical for PDEs than for ODEs (if one does not want to deal with artificial delays represented by delay differential equations). We examined the Kedem-Katchalsky as representing passive transport mechanisms with the difference of concentrations at both sides of the membrane as the driving force for exchanges; however, bigger species are rather transported actively, which should be taken into account in more sophisticated models of nucleocytoplasmic transport [41].
On a more technical note, we have also briefly introduced the semiimplicit Rothe method that can be used for numerical solution of reaction-diffusion systems as an alternative to other used methods; they can be used also in 3D simulations, [31], where other methods are in general more demanding in time and memory.
We hope that the simple presentation of this alternate solution to classical ODEs will be of some help to biologists and modellers who want to describe intracellular spatio-temporal dynamics of proteins in a faithful, yet more demanding in terms of parameter estimation, way. We would as well like to mention that such reaction-diffusion PDE models are also amenable to describe spatio-temporal dynamics at the level of cell populations and that, by introducing intercellular signalling, it is in principle possible to connect the two observation levels. This perspective still remains a challenge to mathematicians and modellers in biology.