Loading [MathJax]/jax/output/SVG/jax.js
Review Recurring Topics

Cognitive Neuroscience of Attention
From brain mechanisms to individual differences in efficiency

  • Aspects of activation, selection and control have been related to attention from early to more recent theoretical models. In this review paper, we present information about different levels of analysis of all three aspects involved in this central function of cognition. Studies in the field of Cognitive Psychology have provided information about the cognitive operations associated with each function as well as experimental tasks to measure them. Using these methods, neuroimaging studies have revealed the circuitry and chronometry of brain reactions while individuals perform marker tasks, aside from neuromodulators involved in each network. Information on the anatomy and circuitry of attention is key to research approaching the neural mechanisms involved in individual differences in efficiency, and how they relate to maturational and genetic/environmental influences. Also, understanding the neural mechanisms related to attention networks provides a way to examine the impact of interventions designed to improve attention skills. In the last section of the paper, we emphasize the importance of the neuroscience approach in order to connect cognition and behavior to underpinning biological and molecular mechanisms providing a framework that is informative to many central aspects of cognition, such as development, psychopathology and intervention.

    Citation: M. Rosario Rueda, Joan P. Pozuelos, Lina M. Cómbita, Lina M. Cómbita. Cognitive Neuroscience of Attention From brain mechanisms to individual differences in efficiency[J]. AIMS Neuroscience, 2015, 2(4): 183-202. doi: 10.3934/Neuroscience.2015.4.183

    Related Papers:

    [1] H.Thomas Banks, Shuhua Hu . Nonlinear stochastic Markov processes and modeling uncertainty in populations. Mathematical Biosciences and Engineering, 2012, 9(1): 1-25. doi: 10.3934/mbe.2012.9.1
    [2] Ying He, Yuting Wei, Junlong Tao, Bo Bi . Stationary distribution and probability density function analysis of a stochastic Microcystins degradation model with distributed delay. Mathematical Biosciences and Engineering, 2024, 21(1): 602-626. doi: 10.3934/mbe.2024026
    [3] Aniello Buonocore, Luigia Caputo, Enrica Pirozzi, Amelia G. Nobile . A non-autonomous stochastic predator-prey model. Mathematical Biosciences and Engineering, 2014, 11(2): 167-188. doi: 10.3934/mbe.2014.11.167
    [4] Mainul Haque, John R. King, Simon Preston, Matthew Loose, David de Pomerai . Mathematical modelling of a microRNA-regulated gene network in Caenorhabditis elegans. Mathematical Biosciences and Engineering, 2020, 17(4): 2881-2904. doi: 10.3934/mbe.2020162
    [5] Dongchen Lu, Wei Han, Kai Lu . Identification of key microRNAs involved in tumorigenesis and prognostic microRNAs in breast cancer. Mathematical Biosciences and Engineering, 2020, 17(4): 2923-2935. doi: 10.3934/mbe.2020164
    [6] Shengjue Xiao, Yufei Zhou, Ailin Liu, Qi Wu, Yue Hu, Jie Liu, Hong Zhu, Ting Yin, Defeng Pan . Uncovering potential novel biomarkers and immune infiltration characteristics in persistent atrial fibrillation using integrated bioinformatics analysis. Mathematical Biosciences and Engineering, 2021, 18(4): 4696-4712. doi: 10.3934/mbe.2021238
    [7] Chris Cosner, Nancy Rodríguez . On the Allee effect and directed movement on the whole space. Mathematical Biosciences and Engineering, 2023, 20(5): 8010-8030. doi: 10.3934/mbe.2023347
    [8] Fan Zhang, Hongtao Liu . Identification of ferroptosis-associated genes exhibiting altered expression in pulmonary arterial hypertension. Mathematical Biosciences and Engineering, 2021, 18(6): 7619-7630. doi: 10.3934/mbe.2021377
    [9] Changxiang Huan, Jiaxin Gao . Insight into the potential pathogenesis of human osteoarthritis via single-cell RNA sequencing data on osteoblasts. Mathematical Biosciences and Engineering, 2022, 19(6): 6344-6361. doi: 10.3934/mbe.2022297
    [10] Guilherme Giovanini, Alan U. Sabino, Luciana R. C. Barros, Alexandre F. Ramos . A comparative analysis of noise properties of stochastic binary models for a self-repressing and for an externally regulating gene. Mathematical Biosciences and Engineering, 2020, 17(5): 5477-5503. doi: 10.3934/mbe.2020295
  • Aspects of activation, selection and control have been related to attention from early to more recent theoretical models. In this review paper, we present information about different levels of analysis of all three aspects involved in this central function of cognition. Studies in the field of Cognitive Psychology have provided information about the cognitive operations associated with each function as well as experimental tasks to measure them. Using these methods, neuroimaging studies have revealed the circuitry and chronometry of brain reactions while individuals perform marker tasks, aside from neuromodulators involved in each network. Information on the anatomy and circuitry of attention is key to research approaching the neural mechanisms involved in individual differences in efficiency, and how they relate to maturational and genetic/environmental influences. Also, understanding the neural mechanisms related to attention networks provides a way to examine the impact of interventions designed to improve attention skills. In the last section of the paper, we emphasize the importance of the neuroscience approach in order to connect cognition and behavior to underpinning biological and molecular mechanisms providing a framework that is informative to many central aspects of cognition, such as development, psychopathology and intervention.



    This paper is concerned with a mathematical model for a gene regulatory network involved in the regulation of DNA transcription. DNA transcription is part of the mechanism by which a sequence of the nuclear DNA is translated into the corresponding protein. The transcription is initiated by the binding of a transcription factor, which is usually another protein, onto the gene's DNA-binding domain. Once bound, the transcription factor promotes the transcription of the nuclear DNA into a messenger RNA (further denoted by mRNA), which, once released, is translated into the corresponding protein by the ribosomes. This process is subject to a high level of noise due to the large variability of the conditions that prevail in the cell and the nucleus at the moment of the transcription. Yet, a rather stable amount of the final protein is needed for the good operation of the cell. The processes that regulate noise levels and maintain cell homeostasis have been scrutinized for a long time. Recently, micro RNAs (further referred to as $ \mu $RNAs) have occupied the front of the scene. These are very short RNAs which do not code for proteins. Many different sorts of $ \mu $RNAs are involved in various epigenetic processes. But one of their roles seems precisely the reduction of noise level in DNA transcription. In this scenario, the $ \mu $RNAs are synthesized together with the mRNAs. Then, some of the synthesized $ \mu $RNAs bind to the mRNAs and de-activate them. These $ \mu $RNA-bound mRNA become unavailable for protein synthesis. It has been proposed that this paradoxical mechanism which seems to reduce the efficiency of DNA transcription may indeed have a role in noise regulation (see [1,2,3] and the review [4]). The goal of the present contribution is to propose a mathematical model of the $ \mu $RNA-mRNA interaction and to use it to investigate the role of $ \mu $RNAs as potential noise regulators.

    Specifically, in this paper, we propose a stochastic chemical kinetic model for the mRNA and $ \mu $RNA content in a cell. The production of mRNAs by the transcription factor and their inactivation through $ \mu $RNA binding are taken into account. More precisely, our model is a simplified version of the circuit used in [5,Figure 2A,A']. We consider a ligand involved in the production of both an mRNA and a $ \mu $RNA, the $ \mu $RNA having the possibility to bind to the mRNA and deactivate it. By contrast to [5], we disregard the way the ligand is produced and consider that the ligand is such that there is a constant production rate of both mRNA and $ \mu $RNA. A second difference to [5] is that we disregard the transcription step of the mRNA into proteins. While [5] proposes to model the $ \mu $RNA as acting on translation, we assume that the $ \mu $RNA directly influences the number of mRNA available for transcription. Therefore, we directly relate the gene expression level to the number of $ \mu $RNA-free mRNA also referred to as the number of unbound mRNA. In order to model the stochastic variability in the production of the RNAs, a multiplicative noise is added to the production rate at all time. From the resulting system of stochastic differential equations, we introduce the joint probability density for mRNA and $ \mu $RNA which solves a deterministic Fokker-Planck equation. The mathematical object of interest is the stationary density solving the Fokker-Planck equation and more precisely the marginal density of the mRNA. The coefficient of variation (also called cell-to-cell variation) of this mRNA density, which is its standard deviation divided by its the expectation, is often considered as the relevant criterion for measuring the robustness of gene expression (see for instance [5]).

    Figure 1.  Exploration of the parameter space. Relative cell to cell variation $ \mathrm{CV}(\rho_\text{fast}^{\delta, \gamma, p}) / \mathrm{CV}(\rho_0^{\delta}) $ for various parameters $ p $, $ \gamma $ and $ \delta $. On the horizontal axis, left means more production of mRNA and right means more production of $ \mu $RNA; On the vertical axis, top means more destruction of mRNA by binding and bottom means more destruction/consumption of mRNA by other mechanisms.
    Figure 2.  Marginal distributions of mRNAs $ \rho_\text{fast}^{\delta, \gamma, p} $ for fast $ \mu $RNAs compared to the free mRNA distribution $ \rho_0^{\delta} $ (black solid curve) for different parameters $ p $ and $ \gamma $. Left: $ \delta = 2 $, $ p = 1.5 $ and $ \gamma $ varies. Right: $ \delta = 2 $, $ \gamma = 1 $ and $ p $ varies.

    Our main goal in this contribution is to provide theoretical and numerical evidence that the robustness of the gene expression is increased in the presence of $ \mu $RNA. At the theoretical level we derive a number of analytical formulas either for particular subsets of parameters of the model or under some time-scale separation hypotheses. From these formulas we can easily compute the cell to cell variation numerically and verify the increased robustness of gene expression when binding with $ \mu $RNA happens in the model. For general sets of parameters, the solution cannot be computed analytically. However we can prove well-posedness of the model and solve the PDE with a specifically designed numerical scheme. From the approximate solution, we compute the coefficient of variation and verify the hypothesis of increased gene expression.

    Another classical approach to the study of noise in gene regulatory networks is through the chemical master equation [6] which is solved numerically by means of Gillespie’s algorithm [7], see e.g., [5,8]. Here, we use a stochastic chemical kinetic model through its associated Kolmogorov-Fokker-Planck equation. Chemical kinetics is a good approximation of the chemical master equation when the number of copies of each molecule is large. This is not the case in a cell where sometimes as few as a 100 copies of some molecules are available. Specifically, including a stochastic term in the chemical kinetic approach is a way to retain some of the randomness of the process while keeping the model complexity tractable. This ultimately leads to a Fokker-Planck model for the joint distribution of mRNAs and $ \mu $RNAs. In [9], a similar chemical kinetic model is introduced with a different modelling of stochasticity. The effect of the noise is taken into account by adding some uncertainty in the (steady) source term and the initial data. The authors are interested in looking at how this uncertainty propagates to the mRNA content and in comparing this uncertainty between situations including $ \mu $RNA production or not. The uncertainty is modeled by random variables with given probability density functions. Compared to [9], the Fokker-Planck approach has the advantage that the random perturbations do not only affect the initial condition and the source term, but are present at all times and vary through time. We believe that this is coherent with how stochasticity in a cell arises through time-varying ecological or biological factors.

    While Fokker-Planck equations are widely used models in mathematical biology [10], their use for the study of gene regulatory network is, up to our knowledge, scarce (see e.g., [11]). Compared to other approaches, the Fokker-Planck model enables us to derive analytical formulas for solutions in certain cases. This is particularly handy for understanding the role of each parameter in the model, calibrating them from real-world data and perform fast numerical computations. Nevertheless, in the general case, the theoretical study and the numerical simulation of the model remains challenging because of the unboundedness of the drift and diffusion coefficients. We believe that we give below all the tools for handling these difficulties, and that our simple model provides a convincing mathematical interpretation of the increase of gene expression in the presence of $ \mu $RNAs.

    The paper is organized as follows. In Section 2, we introduce the system of SDEs and the corresponding Fokker-Planck models. In Section 3, we discuss the well-posedness of the Fokker-Planck equations and derive analytical formulas for solutions under some simplifying hypotheses. In Section 4, we use the analytical formulas for solutions to give mathematical and numerical proofs of the decrease of cell-to-cell variation in the presence of $ \mu $RNA. In Section 5, we propose a numerical scheme for solving the main Fokker-Planck model and gather further evidence confirming the hypothesis of increased gene expression from the simulations. Finally, in Section 6 we discuss the particular choice of multiplicative noise (i.e., the diffusion coefficient in the Fokker-Planck equation) in our model. In the appendix, we derive weighted Poincaré inequalities for gamma and inverse-gamma distributions which are useful in the analysis of Section 3. The code used for numerical simulations in this paper is publicly available on GitLab [12].

    In this section, we introduce three steady Fokker-Planck models whose solutions describe the distribution of unbound mRNA and $ \mu $RNA within a cell. The solutions to these equations can be interpreted as the probability density functions associated with the steady states of stochastic chemical kinetic systems describing the production and destruction of mRNA and $ \mu $RNA. In Section 2.1 we introduce the main model for which the consumption of RNAs is either due to external factors in the cell (translation, {etc.}) or to binding between the two types of mRNA and $ \mu $RNA. Then, for comparison, in Section 2.2 we introduce the same model without binding between RNAs. Finally in Section 2.3, we derive an approximate version of the first model, by considering that reactions involving $ \mu $RNAs are infinitely faster than those involving mRNAs, which amplifies the binding phenomenon and mathematically allows for the derivation of analytical formulas for solutions. The latter will be made explicit in Section 3.

    We denote by $ r_t $ the number of unbound mRNA and $ \mu_t $ the number of unbound $ \mu $RNA of a given cell at time t. The kinetics of unbound mRNA and $ \mu $RNA is then given by the following stochastic differential equations

    $ {drt=(crcrtμtkrrt)dt+ 2σrrtdB1t,dμt=(cμcrtμtkμμt)dt+2σμμtdB2t, $ (2.1)

    with $ c_r $, $ c_\mu $, $ k_r $, $ k_\mu $, $ \sigma_r $, $ \sigma_\mu $ being some given positive constants and $ c $ being a given non-negative constant. Let us detail the meaning of each term in the modeling. The first term of each equation models the constant production of mRNA (resp. $ \mu $RNA) by the ligand at a rate $ c_r $ (resp. $ c_\mu $). The second term models the binding of the $ \mu $RNA to the mRNA. Unbound mRNA and $ \mu $RNA are consumed by this process at the same rate. The rate increases with both the number of mRNA and $ \mu $RNA. In the third term, the parameters $ k_r $ and $ k_\mu $ are the rates of consumption of the unbound mRNA or $ \mu $RNA by various decay mechanisms. The last term in both equations represents stochastic fluctuations in the production and destruction mechanisms of each species. It relies on a white noise $ \, \mathrm{d} B_t / \, \mathrm{d} t $ where $ B_t = (B_t^1, B_t^2) $ is a two-dimensional standard Brownian motion. The intensity of the stochastic noise is quantified by the parameters $ \sqrt{2\sigma_r}\, r_t $ and $ \sqrt{2\sigma_\mu}\, \mu_t $. Such a choice of multiplicative noise ensures that $ r_t $ and $ \mu_t $ remain non-negative along the dynamics. The Brownian motions $ B_t^1 $ and $ B_t^2 $ are uncorrelated. The study of correlated noises or the introduction of extrinsic noise sources would be interesting, but will be discarded here.

    In this paper we are interested in the invariant measure of Eq (2.1) rather than the time dynamics described by the above SDEs. From the modelling point of view, we are considering a large number of identical cells and we assume that mRNA and $ \mu $RNA numbers evolve according to Eq (2.1). Then we measure the distribution of both RNAs among the population, when it has reached a steady state $ f \equiv f(r, \mu) $. According to Itô's formula, the steady state should satisfy the following steady Fokker-Planck equation

    $ {Lf(r,μ) = 0,(r,μ)Ω=(0,)2,Ωf(r,μ)drdμ = 1,f(r,μ)  0. $ (2.2)

    where the Fokker-Planck operator is given by

    $ Lf(r,μ) := r[r(σrr2f)(crcrμkrr)f]+μ[μ(σμμ2f)(cμcrμkμμ)f]. $ (2.3)

    Since we do not model the protein production stage, we assume that the observed distribution of gene expression level is proportional to the marginal distribution of mRNA, i.e..

    $ \rho(r) = \int_0^\infty f(r, \mu) \, \mathrm{d} \mu\, . $

    By integration of Eq (2.2) in the $ \mu $ variable, $ \rho $ satisfies the equation

    $ r[r(σr r2 ρ)( crc r jμ(r)kr r)ρ]=0. $ (2.4)

    The quantity $ j_\mu(r) $ is the conditional expectation of the number of $ \mu $RNA within the population in the presence of $ r $ molecules of mRNA and it is given by

    $ jμ(r)=1ρ(r)μ f(r,μ)dμ. $ (2.5)

    Before ending this paragraph, we note an alternate way to derive the Fokker-Planck Eq (2.2) from the chemical master equation through the chemical Langevin equation. We refer the interested reader to [13].

    In the case where there is no $ \mu $RNA binding, namely when $ c = 0 $, the variables $ r $ and $ \mu $ are independent. Thus, the densitites of the invariant measures satisfying Eq (2.2) are of the form

    $ f_0(r, \mu) = \rho_0(r)\lambda(\mu)\, , $

    where $ \lambda(\mu) $ is the density of the marginal distrubution of $ \mu $RNA. From the modelling point of view, it corresponds to the case where there is no feed-forward loop from $ \mu $RNA. Therefore, only the dynamics on mRNA, and thus $ \rho_0 $, is of interest in our study. It satisfies the following steady Fokker-Planck equation obtained directly from Eq (2.4),

    $ {r[σrr(r2ρ0)(crkrr)ρ0] = 0,0ρ0(r)dr = 1,ρ0(r)0. $ (2.6)

    It can be solved explicitly as we will discuss in Section 3.2.

    The Fokker-Planck Eq (2.2) cannot be solved explicitly. However, one can make some additional assumptions in order to get an explicit invariant measure providing some insight into the influence of the binding mechanism with $ \mu $RNA. This is the purpose of the model considered hereafter.

    Let us assume the $ \mu $RNA-mRNA binding rate, the $ \mu $RNA decay and the noise on $ \mu $RNA are large. Since the sink term of the $ \mu $RNA equation is large, it is also natural to assume that the $ \mu $RNA content is small. Mathematically, we assume the following scaling

    $ c = \frac{\tilde{c}}{\varepsilon}\, , \quad k_\mu = \frac{\tilde{k}_\mu}{\varepsilon}\, , \quad \sigma_\mu = \frac{\tilde{\sigma}_\mu}{\varepsilon}\, , \quad \mu_t = \varepsilon\tilde{\mu_t}\, , $

    for some small constant $ \varepsilon > 0 $. Then $ (r_t, \tilde{\mu}_t) $ satisfies

    $ \left\{ drt=(cr˜crt˜μt˜krrt)dt+ 2σrrtdB1t,εd˜μt=(cμ˜crt˜μt˜kμ˜μt)dt+2ε˜σμ˜μtdB2t, \right. $

    whose corresponding steady Fokker-Planck equation for the invariant measure then writes, dropping the tilde,

    $ \partial_r\left[\partial_r (\sigma_rr^2f_\varepsilon) -(c_r - c\, r\, \mu - k_r\, r)f_\varepsilon\right]+\frac{1}{\varepsilon} \partial_\mu\left[ \partial_\mu( \sigma_\mu\mu^2f_\varepsilon) -(c_\mu - c\, r\, \mu - k_\mu\, \mu)f_\varepsilon\right]\ = \ 0. $

    In the limit case where $ \varepsilon\to 0 $, one may expect that at least formally, the density $ f_\varepsilon $ converges to a limit density $ f_\text{fast} $ satisfying

    $ \partial_\mu\left[ \partial_\mu( \sigma_\mu\mu^2f_\text{fast}) -(c_\mu - c\, r\, \mu - k_\mu\, \mu)f_\text{fast}\right]\ = \ 0\, . $

    As $ r $ is only a parameter in the previous equation and since the first marginal of $ f_\varepsilon $ still satisfies Eq (2.4) for all $ \varepsilon $, one should have (formally)

    $ {ffast(r,μ) = ρfast(r)M(r,μ)  0,μ[μ(σμμ2M)(cμcrμkμμ)M] = 0,r[r(σr r2 ρfast)( crc r jfast(r)kr r)ρfast]=0,0ρfast(r)dr = 1, 0M(r,μ)dμ = 1,jfast(r) = 0μ M(r,μ)dμ. $ (2.7)

    In this section, we show that the three previous models are well-posed. For the Fokker-Planck Eqs (2.6) and (2.7), we explicitly compute the solutions. They involve inverse gamma distributions.

    The expressions of the gamma and inverse gamma probability densities are respectively

    $ γα,β(x) = Cα,βxα1exp(βx), $ (3.1)

    and

    $ gα,β(y) = Cα,βy1+αexp(βy), $ (3.2)

    for $ x, y\in(0, \infty) $. The normalization constant is given by $ C_{\alpha, \beta} = \beta^{\alpha}/\Gamma(\alpha) $ where $ \Gamma $ is the Gamma function. Observe that by the change of variable $ y = 1/x $ one has

    $ g_{\alpha, \beta}(y)\, \mathrm{d} y\ = \ \gamma_{\alpha, \beta}(x)\, \mathrm{d} x $

    which justifies the terminology. Let us also recall that the first and second moments of the inverse gamma distribution are

    $ 0ygα,β(y)dy = βα1, if α>1, β>0, $ (3.3)
    $ 0y2gα,β(y)dy = β2(α1)(α2), if α>2, β>0, $ (3.4)

    Interestingly enough, we can show (see Appendix A.1 for details and additional results) that inverse gamma distributions with finite first moment ($ \alpha > 1 $) satisfy a (weighted) Poincaré inequality. The proof of the following proposition is done in Appendix A.1 among more general considerations.

    Proposition 3.1. Let $ \alpha > 1 $ and $ \beta > 0 $. Then, for any function $ v $ such that the integrals make sense, one has

    $ 0|v(y)vgα,β|2gα,β(y)dy  1α10|v(y)|2gα,β(y)y2dy, $ (3.5)

    where for any probability density $ \nu $ and any function $ u $ on $ (0, \infty) $, the notation $ \langle\, u\, \rangle_{\nu} $ denotes $ \int u\nu $.

    In the case of free mRNAs, a solution to Eq (2.6) can be computed explicitely and takes the form of an inverse gamma distribution.

    Lemma 3.2. The following inverse gamma distribution

    $ ρ0(r) = g1+krσr,crσr(r) = C1+krσr,crσr1r2+krσrexp(crσrr) $ (3.6)

    is the only classical solution to Eq (2.6).

    Proof. First observe that

    $ \partial_r\left[\sigma_r\partial_r (r^2\rho_0) -(c_r - k_r\, r)\rho_0\right]\ = \ \partial_r\left[\sigma_r\, r^2\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}\partial_r \left(\rho_0\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1}\right)\right]\, . $

    Therefore a solution of Eq (3.1) must be of the form

    $ \rho_0(r)\ = \ C_1\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}\int_{1}^r\sigma_r^{-1}r^{-2}g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1}(r)\, \mathrm{d} r + C_2\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}\, , $

    for some constants $ C_1, C_2 $. The first term decays like $ 1/r $ at infinity, thus the only probability density $ \rho_0 $ of this form is obtained for $ C_1 = 0 $ and $ C_2 = 1 $.

    The Poincaré inequality (3.5) tells us that the solution of Lemma 3.2 is also the only (variational) solution in the appropriate weighted Sobolev space. Indeed, we may introduce the natural Hilbert space associated with Eq (2.6),

    $ X_{\alpha, \beta}\ = \ \{v:(0, \infty)\to\mathbb{R}\, , \ \|v\|_{X_{\alpha, \beta}} \lt \infty\} $

    with a squared norm given by

    $ \|v\|_{X_{\alpha, \beta}}^2\ = \ \int_0^\infty \left(\left[\left(v/g_{\alpha, \beta}\right)(y)\right]^2 + \left[(v/g_{\alpha, \beta})'(y)\right]^2\, y^2\right)\, g_{\alpha, \beta}(y)\, \mathrm{d} y\, . $

    Then the following uniqueness result holds.

    Lemma 3.3. The classical solution $ \rho_0 = g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}} $ is the only solution of Eq (2.6) in $ X_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}} $.

    Proof. If $ \rho_0 $ and $ \tilde{\rho_0} $ are two solutions of Eq (2.6), a straightforward consequence of inequality (3.5) is that $ \|\rho_0-\tilde{\rho_0}\|_{X_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}} = 0 $. This is obtained by integrating the difference between the equation on $ \rho_0 $ and $ \tilde{\rho_0} $ against $ (\rho_0-\tilde{\rho_0})g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1} $.

    Another consequence of the Poincaré inequality is that if we consider the time evolution associated with Eq (2.6) then solutions converge exponentially fast towards the steady state $ \rho_0 $. This justifies our focus on the stationary equations. The transient regime is very short and equilibrium is reached quickly. We can quantify the rate of convergence in terms of the parameters.

    Proposition 3.4. Let $ \xi $ solve the Fokker-Planck equation

    $ \partial_t\xi(t, r)\ = \ \partial_r\left[\sigma_r\partial_r (r^2\xi(t, r)) -(c_r - k_r\, r)\xi(t, r)\right]\, , $

    starting from the probability density $ \xi(0, r, \mu) = \xi^\mathit{\text{in}}(r, \mu) $. Then for all $ t\geq0 $,

    $ \int_0^\infty(\xi(t, r)-\rho_0(r))^2\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1}(r)\, \mathrm{d} r \leq\ e^{-\frac{k_r}{\sigma_r}\, t}\, \int_0^\infty(\xi^\mathit{\text{in}}(r)-\rho_0(r))^2\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1}(r)\, \mathrm{d} r\, . $

    Proof. Observe that $ \xi-\rho_0 $ solves the unsteady Fokker-Planck equation, so that by multiplying the equation by $ (\xi-\rho_0)\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1} $ and integrating in $ r $ one gets

    $ \frac{\, \mathrm{d}}{\, \mathrm{d} t} \int_0^\infty(\xi(t, r)-\rho_0(r))^2\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}^{-1}(r)\, \mathrm{d} r \, +\, \int_0^\infty\left|\partial_r\left(\frac{\xi(t, \cdot)-\rho_0}{g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}}\right)(r)\right|^2\, g_{1+\frac{k_r}{\sigma_r}, \frac{c_r}{\sigma_r}}(r)\, r^2\, \, \mathrm{d} r\ = \ 0\, . $

    Then by using the Poincaré inequality Eq (3.5) and a Gronwall type argument, one gets the result.

    Now we focus on the solution of Eq (2.7). The same arguments as those establishing Lemma 3.3 show that the only function $ M $ satisfying Eq (2.7) is the following inverse gamma distribution

    $ M(r,μ) = g1+kμσμ+cσμr,cμσμ(μ). $ (3.7)

    Then an application of Eq (3.3) yields

    $ jfast(r)=cμkμ+cr. $ (3.8)

    It remains to find $ \rho_\text{fast} $ which is a probability density solving the Fokker-Planck equation

    $ \partial_r\left[\partial_r(\sigma_r\ r^2\ \rho_\text{fast}) - (\ c_r - \frac{c_\mu\, c\, r}{k_\mu+cr} - k_r\ r)\rho_\text{fast} \right] = 0\, . $

    Arguing as in the proof of Lemma 3.2, one observe that integrability properties force $ \rho_\text{fast} $ to actually solve

    $ \partial_r(\sigma_r\ r^2\ \rho_\text{fast}) - (\ c_r - \frac{c_\mu\, c\, r}{k_\mu+cr} - k_r\ r)\rho_\text{fast} = 0\, . $

    which yields

    $ ρfast(r) = C(1+kμcr)ccμσrkμ1r2+krσrexp(crσrr) $ (3.9)

    where $ C\equiv C(c_r, c_\mu, k_r, k_\mu, \sigma_r, \sigma_\mu, c) $ is a normalizing constant making $ \rho_\text{fast} $ a probability density function.

    Now we are interested in the well-posedness of Eq (2.2), for which we cannot derive explicit formulas anymore. Despite the convenient functional framework introduced in Section 3.2, classical arguments from elliptic partial differential equation theory do not seem to be adaptable to the case $ c > 0 $. The main obstruction comes from an incompatibility between the natural decay of functions in the space $ X_{\alpha, \beta} $ and the rapid growth of the term $ c\, r\, \mu $ when $ |(r, \mu)|\to \infty $.

    However, thanks to the results of [14] focused specifically on Fokker-Planck equations, we are able to prove well-posedness of the steady Fokker-Planck Eq (2.2). The method is based on finding a Lyapunov function for the adjoint of the Fokker-Planck operator and relies on an integral identity proved by the same authors in [15]. The interested reader may also find additional material and a comprehensive exposition concerning the analysis of general Fokker-Planck equations for measures in [16].

    First of all let us specify the notion of solution. A weak solution to Eq (2.2) is an integrable function $ f $ such that

    $ {Ωf(r,μ)Lφ(r,μ) = 0,for all φCc(Ω),Ωf(r,μ) = 1,f(r,μ)  0, $ (3.10)

    where the adjoint operator is given by

    $ Lφ(r,μ) := σrr22rrφ+(crcrμkrr)rφ+σμμ22μμφ+(cμcrμkμμ)μφ. $ (3.11)

    A reformulation and combination of [14,Theorem A and Proposition 2.1] provides the following result.

    Proposition 3.5 ([14]). Assume that there is a smooth function $ U:\Omega\rightarrow[0, +\infty) $, called Lyapunov function with respect to $ \mathcal{L} $, such that

    $ lim(r,μ)¯ΩU(r,μ) = +, $ (3.12)

    and

    $ lim(r,μ)¯ΩLU(r,μ) = , $ (3.13)

    where $ \overline{\partial\Omega} = \partial\Omega\cup(\{+\infty\}\times\mathbb{R}_+)\cup(\mathbb{R}_+\times\{+\infty\}) $. Then there is a unique $ f $ satisfying Eq (3.10). Moreover $ f\in W^{1, \infty}_\mathrm{loc}(\Omega) $.

    Remark 3.6. The method of Lyapunov functions is a standard tool for proving well-posedness of many problems in the theory of ordinary differential equations, dynamical systems... For diffusion processes and Fokker-Planck equations its use dates back to Has'minskii [17]. We refer to [16,Chapter 2], [14] and references therein for further comments on the topic. Let us stress however that Lyapunov functions are not related (at least directly) to the Lyapunov (or entropy) method for evolution PDEs in which one shows the monotony of a functional to quantify long-time behavior.

    Remark 3.7. Thanks of the degeneracy of the diffusivities at $ r = 0 $ and $ \mu = 0 $ and the Lyapunov function condition, one doesn't need supplementary boundary conditions in Eq (3.10) for the problem to have a unique solution. This is different from standard elliptic theory where boundary conditions are necessary to define a unique solution when the domain and the coefficients are bounded with uniformly elliptic diffusivities. Further comments may be found in [14].

    Lemma 3.8. Choose any two constants $ b_r > \frac{c}{k_\mu} $ and $ b_\mu > \frac{c}{k_r} $. Then, the function $ U:\Omega\rightarrow\mathbb{R} $ defined by

    $ U(r, \mu)\ = \ b_rr-\ln(b_rr)+b_\mu\mu-\ln(b_\mu\mu) $

    is a Lyapunov function with respect to $ \mathcal{L} $ (i.e., it is positive on $ \Omega $ and it satisfies Eqs (3.12) and (3.13)).

    Proof. First observe that condition (3.12) is clearly satisfied. Also, $ U $ is minimal at $ (b_r^{-1}, b_\mu^{-1}) $ where it takes the value $ 2 $ and thus it is positive on $ \Omega $. Finally a direct computation yields

    $ \mathcal{L}U(r, \mu)\ = \ (\sigma_r+\sigma_\mu+b_rc_r+b_\mu c_\mu+k_r+k_\mu) -\frac{c_r}{r}-\frac{c_\mu}{\mu} -(b_rk_r-c)r-(b_\mu k_\mu-c)\mu-c r\mu(b_r+b_\mu)\, , $

    and Eq (3.13) follows. Now we state our well-posedness result for the Fokker-Planck Eq (3.10).

    Proposition 3.9. There is a unique weak solution $ f $ to the steady Fokker-Planck Eq (2.2). Moreover, $ f $ is indefinitely differentiable in $ \Omega $.

    Proof. The existence and uniqueness of a solution $ f\in W^{1, \infty}_\mathrm{loc}(\Omega) $ is a combination of Proposition 3.5 and Lemma 3.8. From there in any smooth compact subdomain $ K\subset\subset\Omega $, we get from standard elliptic theory [18] that $ f\in\mathcal{C}^\infty(K) $, since the coefficents are smooth and the operator is uniformly elliptic.

    In this section we focus on the comparison between the explicit distributions given by Eqs (3.6) and (3.9). We are providing theoretical and numerical evidence that the coefficient of variation (which is a normalized standard deviation) of distribution (3.9) is less than that of distribution (3.6). This quantity called cell to cell variation in the biological literature [5] characterizes the robustness of the gene expression level (the lower the better). We start by performing a rescaling in order to extract the dimensionless parameters which characterize the distributions.

    In order to identify the parameters of importance in the models, we rescale the variables $ r $ and $ \mu $ around characteristic values $ \bar{r} $ and $ \bar{\mu} $ chosen to be

    $ ˉr=crkrandˉμ = cμkμ. $ (4.1)

    These choices are natural in the sense that they correspond to the steady states of the mRNA and $ \mu $RNA dynamics without binding nor stochastic effects, that is respectively $ \, \mathrm{d} r_t = (c_r - k_r\ r_t)\, \mathrm{d} t $ and $ \, \mathrm{d} \mu_t = (c_\mu - k_\mu\ \mu_t)\, \mathrm{d} t $. When the noise term is added, it still corresponds to the expectation of the invariant distribution, that is the first moment of $ \rho_0 $ in the case of mRNA. We introduce $ f^\text{ad} $ such that for all $ (r, \mu)\in\Omega $ one has

    $ \frac{1}{\bar{r}\bar{\mu}}\, f^\text{ad}\left(\frac{r}{\bar{r}}, \frac{\mu}{\bar{\mu}}\right)\ = \ f(r, \mu)\, . $

    After some computations one obtains that the Fokker-Planck Eqs (2.2) and (2.3) can be rewritten in terms of $ f^\text{ad} $ as

    $ r[δ(1γprμr)fadr(r2fad)]+ μ[δκ(1γrμμ)fadνμ(μ2fad)] = 0, $ (4.2)

    The marginal distributions $ \rho_0 $ and $ \rho_\text{fast} $ are rescaled into dimensionless densities

    $ ρδ0(r) = g1+δ,δ(r) = Cad01r2+δexp(δr) $ (4.3)
    $ ρδ,γ,pfast(r) = Cadfast(1+1γr)γpδ1r2+δexp(δr) $ (4.4)

    where $ C_0^\text{ad} $ and $ C_\text{fast}^\text{ad} $ are normalizing constants depending on the parameters of the model and $ \delta $, $ p $ and $ \gamma $ are dimensionless parameters. The first parameter

    $ δ = krσr, $ (4.5)

    only depends on constants that are independent of the dynamics of $ \mu $RNAs. The two other dimensionless parameters appearing in the marginal distribution of mRNA in the presence of fast $ \mu $RNA are

    $ p = cμcr, $ (4.6)

    and

    $ γ = cˉrkμ = ccrkμkr. $ (4.7)

    Let us give some insight into the biological meaning of these parameters. The parameter $ \gamma $ measures the relative importance of the two mechanisms of destruction of $ \mu $RNAs, namely the binding with mRNAs versus the natural destruction/consumption. A large $ \gamma $ means that the binding effect is strong and conversely. The parameter $ p $ compares the production rate of $ \mu $RNAs with that of mRNAs. Large values of $ p $ mean that there are much more $ \mu $RNAs than mRNAs produced per unit of time.

    Finally, in the Fokker-Planck model (4.2), there are also two other parameters which are

    $ κ = kμkr, $ (4.8)

    and

    $ ν = σμσr. $ (4.9)

    The parameter $ \kappa $ compares consumption of $ \mu $RNA versus that of mRNA by mechanisms which are not the binding between the two RNAs. The parameter $ \nu $ compares the amplitude of the noise in the dynamics of $ \mu $RNA versus that of the mRNA.

    Remark 4.1. Observe that the approximation of fast $ \mu $RNA leading to the model discussed in Section 2.3 in its dimensionless form amounts to taking $ \nu = \kappa = 1/\varepsilon $ and letting $ \varepsilon $ tend to $ 0 $.

    For any suitably integrable non-negative function $ \nu $, let us denote by

    $ m_k(\nu)\ = \ \int y^k\, \nu(y)\, \mathrm{d} y $

    its $ k $-th moment. The coefficient of variation or cell to cell variation (CV) is defined by

    $ CV(ν) = Var(ν/m0(ν))1/2Exp(ν/m0(ν)) = (m2(ν)m0(ν)m1(ν)21)1/2 $ (4.10)

    where $ \mathrm{Exp}(\cdot) $ and $ \mathrm{Var}(\cdot) $ denote the expectation and variance. Let us state a first lemma concerning some cases where the coefficient of variation can be computed exactly.

    Lemma 4.2. Consider the dimensionless distributions defined in Eqs (4.3) and (4.4). Then one has that

    $ Exp(ρδ0) = 1,Var(ρδ0) = 1δ1,CV(ρδ0) = 1δ1, $ (4.11)

    where the variance and coefficient of variation are well-defined only for $ \delta > 1 $. Then for any $ \delta > 1 $, the following limits holds

    $ limγ0CV(ρδ,γ,pfast)=CV(ρδ0),p>0,limp0CV(ρδ,γ,pfast)=CV(ρδ0),γ>0,limγCV(ρδ,γ,pfast)=CV(ρδ0),p[0,1). $

    Proof. The formulas for the moments follow from Eqs (3.3) and (3.4). Then observe that for all $ r $, one has

    $ \lim\limits_{\gamma\to0}\rho_\text{fast}^{\delta, \gamma, p}(r)\ = \ \lim\limits_{p\to0}\rho_\text{fast}^{\delta, \gamma, p}(r)\ = \rho_0^{\delta}(r) $

    and

    $ \lim\limits_{\gamma\to+\infty}\rho_\text{fast}^{\delta, \gamma, p}(r)\ = \ g_{1+\delta, (1-p)\delta}(r) $

    and one can then take limits in integrals by dominated convergence.

    Let us give a biological interpretation of the previous lemma. When $ \gamma = 0 $ or $ p = 0 $, which respectively corresponds to the cases where there is no binding between mRNA and $ \mu $RNA or there is no production of $ \mu $RNA, the coefficient of variation is unchanged from the case of free mRNAs. The last limit states that if the $ \mu $RNA production is weaker than the mRNA production, then in the regime where all $ \mu $RNA is consumed by binding with mRNA, the coefficient of variation is also unchanged.

    Outside of these asymptotic regimes, the theoretical result one would like to have is the following.

    Conjecture 4.3. For any $ \delta > 1 $, $ \gamma, p > 0 $ and one has $ \mathrm{CV}(\rho_\mathit{\text{fast}}^{\delta, \gamma, p})\ \leq\ \mathrm{CV}(\rho_0^{\delta}) $.

    At the moment, we are able to obtain the following uniform in $ \gamma $ and $ p $ bound

    $ CV(ρδ,γ,pfast)  Cδ := ((δδ1)2(11(δ1)2)δ21)12, $ (4.12)

    which holds for all $ \delta > 2 $, $ \gamma > 0 $ and $ p\geq0 $. The result is proved in Proposition in the Appendix. Observe that $ C_\delta\geq \mathrm{CV}(\rho_0^{\delta}) $ but asymptotically

    $ C_\delta\sim_{\delta\to\infty}\mathrm{CV}(\rho_0^{\delta})\ = \ \frac{1}{\sqrt{\delta-1}}\, , $

    so that $ C_\delta $ is fairly close to $ \mathrm{CV}(\rho_0^{\delta}) $ for large $ \delta $.

    In the next section we provide numerical evidence that it should be possible to improve the right-hand side of the bound (4.12) and prove Conjecture 4.3. Let us also mention that using integration by parts formulas it is possible to establish a recurrence relation between moments. From there one can infer the inequality of Conjecture 4.3 for subsets of parameters $ (\gamma, p) $. As the limitation to these subsets is purely technical and do not have any particular biological interpretation we do not report these results here.

    Now, we explore the space of parameters $ (\delta, \gamma, p) $ in order to compare the cell to cell variation in the case of fast $ \mu $RNA and in the case of free mRNA.

    In order to evaluate numerically the cell to cell variation we need to compute $ m_k(\rho_\text{fast}^{\delta, \gamma, p}) $, for $ k = 0, 1, 2 $. Observe that after a change of variable these quantities can be rewritten (up to an explicit multiplicative constant depending on parameters)

    $ I_k = \int_0^\infty f_k(s)\, s^{\delta-2}e^{-s}\, \mathrm{d} s\, , $

    with $ f_k(s) = s^{2-k}(1+s/(\gamma\delta))^{p\gamma\delta} $. For the numerical computation of these integrals, we use a Gauss-Laguerre quadrature

    $ I_k\ \approx\ \sum\limits_{i = 1}^N\omega_i^Nf_k(x_i^N)\, . $

    which is natural and efficient as we are dealing with functions integrated against a gamma distribution. We refer to [19] and references therein for the definition of the coefficients $ \omega_i^N $ and quadrature points $ x_i^N $. The truncation order $ N $ is chosen such that the numerical error between the approximation at order $ N $ and $ N+1 $ is inferior to the given precision $ 10^{-8} $ when $ p\leq1 $. For $ p\geq1 $, the function $ f_k $ may take large values and it is harder to get the same numerical precision. In the numerical results below the mean error for the chosen sets of parameters with large values of $ p $ is around $ 10^{-4} $ and the maximal error is $ 10^{-2} $. This is good enough to comment on qualitative behavior. The code used for these numerical simulations is publicly available on GitLab [12].

    We plot the relative cell to cell variation $ \mathrm{CV}(\rho_\text{fast}^{\delta, \gamma, p}) / \mathrm{CV}(\rho_0^{\delta}) $ with respect to $ \gamma $ and $ p $ for two different values of $ \delta $. The results are displayed on Figure 1. Then, on Figure 2, we draw the explicit distributions $ \rho_\text{fast}^{\delta, \gamma, p} $ for various sets of parameters and compare it with $ \rho_0^{\delta} $.

    The numerical simulations of Figure 1 suggest that the bound (4.12) is non-optimal and Conjecture 4.3 should be satisfied. Observe also that the asymptotics of Lemma 4.2 are illustrated.

    From a modeling point of view, these simulations confirm that for any choice of parameter, the presence of (fast) $ \mu $RNA makes the cell to cell variation decrease compared to the case without $ \mu $RNA. Moreover, the qualitative behavior with respect to the parameters makes sense. Indeed we observe that whenever enough $ \mu $RNA is produced ($ p\geq1 $), the increase of the binding phenomenon ($ \gamma\to\infty $) makes the cell to cell variation decay drastically.

    In this section, we compute the gene expression level of the main model described by Eq (2.2). In this case, as there is no explicit formula for the solution, we will compute an approximation of it using a discretization of the Fokker-Planck equation. In order to compute the solution in practice, we restrict the domain to the bounded domain $ \Omega_b = [r_{\min}, r_{\max}]\times[\mu_{\min}, \mu_{\max}] $. Because of the truncation, we add zero-flux boundary conditions in order to keep a conservative equation. It leads to the problem

    $ {r[(crcrμkrr)fr(σrr2f)]+ μ[(cμcrμkμμ)fμ(σμμ2f)] = 0,in Ωbr(σrr2f)(crcrμkrr)f=0, if r=rmin or rmax,μ(σμμ2f)(cμcrμkμμ)f=0, if μ=μmin or μmax,Ωbfdrdμ = 1.  $ (5.1)

    In order for the numerical scheme to be more robust with respect to the size of the parameters, we discretize the equation in dimensionless version according to Eq (4.2). It will also allow for comparisons with numerical experiments of the previous sections.

    As the coefficients in the advection and diffusion parts of Eq (4.2) grow rapidly in $ r $, $ \mu $ and degenerate when $ r = 0 $ and $ \mu = 0 $, the design of an efficient numerical solver for Eq (4.2) is not straightforward. Moreover a desirable feature of the scheme would be a preservation of the analytically known solution corresponding to $ \gamma = 0 $. Because of these considerations we will discretize a reformulated version of the equation in which the underlying inverse gamma distributions explicitly appear. It will allow for a better numerical approximation when $ r $ and $ \mu $ are either close to $ 0 $ or large. The reformulation is the following

    $ r[r2h(1)(r,μ)r(fadh(1)(r,μ))]μ[νμ2h(2)(r,μ)μ(fadh(2)(r,μ))] = 0, $ (5.2)

    with the associated no-flux boundary conditions and where the functions $ h^{(1)} $ and $ h^{(2)} $ are given by

    $ h(1)(r,μ) = r(1+pμγ)δ2exp(δr), $ (5.3)

    and

    $ h(2)(r,μ) = μ(1+rγ)δκν2exp(δκνr). $ (5.4)

    We use a discretization based on the reformulation (5.2). It is inspired by [20] and is fairly close to the so-called Chang-Cooper scheme [21].

    We use a finite-volume scheme. The rectangle $ \Omega_b $ is discretized with a structured regular mesh of size $ \Delta r $ and $ \Delta \mu $ in each respective direction. The centers of the control volumes are the points $ \left(r_i, \mu_j\right) $ with $ r_i = \Delta r/2 + i\Delta r $ and $ \mu_j = \Delta \mu/2 + j\Delta \mu $ for $ i\in\{0, \dots, N_r-1\} $ and $ j\in\{0, \dots, N_\mu-1\} $. We also introduce the intermediate points $ r_{i+1/2} $ with $ i\in\{-1, \dots, N_r-1\} $ and $ \mu_{j+1/2} $ with $ j\in\{-1, \dots, N_\mu-1\} $ defined with the same formula as before. The approximation of the solution on the cell $ (i, j) $ is denoted by

    $ f_{ij} \approx \frac{1}{\Delta r \Delta \mu} \int_{r_{i-1/2}}^{r_{i+1/2}}\int_{\mu_{j-1/2}}^{\mu_{j+1/2}} f^\text{ad}(r, \mu) \, \mathrm{d} r \, \mathrm{d} \mu. $

    The scheme reads, for all $ i\in\{0, \dots, N_r-1\} $ and $ j\in\{0, \dots, N_\mu-1\} $,

    $ {Fi+1/2,jFi1/2,j+Gi,j+1/2Gi,j1/2=0,FNr1/2,j = F1/2,j = Gi,Nμ1/2 = Gi,1/2 = 0i,jfijΔrΔμ = 1 $ (5.5)

    where the fluxes are given by a centered discretization of the reformulation (5.2), namely

    $ Fi+1/2,j=ΔμΔrr2i+1/2(h(1)(ri+1/2,μj)h(1)(ri+1,μj)fi+1,jh(1)(ri+1/2,μj)h(1)(ri,μj)fij), $ (5.6)

    and

    $ Gi,j+1/2=νΔrΔμμ2j+1/2(h(2)(ri,μj+1/2)h(2)(ri,μj+1)fi,j+1h(2)(ri,μj+1/2)h(2)(ri,μj)fij). $ (5.7)

    One can show that the scheme (5.5) possesses a unique solution which is non-negative by following, for instance, the arguments of [22,Proposition 3.1]. Moreover, by construction, the scheme is exact in the case $ \gamma = 0 $.

    Remark 5.1 (Choice of $ r_{\min} $, $ r_{\max} $, $ \mu_{\min} $, $ \mu_{\max} $). Clearly $ f $ decays faster at infinity than $ \rho_0 $ since the convection term coming from the binding phenomenon brings mass closer to the origin. Therefore an appropriate choice for $ r_{\max} $ and $ \mu_{\max} $, coming from the decay of the involved inverse gamma distributions, should be (say) $ r_{\max}^{-\delta}\leq 10^{-8} $ and $ \mu_{\max}^{-\delta\kappa/\nu}\leq 10^{-8} $ so that the error coming from the tails of the distributions in the computation of moments is negligible. Similarly, near the origin the distributions decay very quickly to $ 0 $ (as $ \exp(-1/\cdot) $). Therefore $ \mu_{\min} $, $ r_{\min} $ can be taken not too small without influencing the precision in the computation of moments of the solution. In practice, we chose $ \mu_{\min} = r_{\min} = 0.06 $. Observe that even if nothing prevents the choice $ \mu_{\min} = r_{\min} = 0 $ on paper, one experiences in practice a bad conditioning of the matrix which has to be inverted for solving the scheme.

    Remark 5.2 (Implementation). Observe that the matrix which has to be inverted in order to solve the scheme is not a square matrix because of the mass constraint (which is necessary to ensure uniqueness of the solution). In practice, in order to solve the corresponding linear system $ MF = B $ where $ F = (f_{ij})_{ij} $ and $ B = (0, \dots, 0, 1)\in\mathbb{R}^{N_rN_\mu+1} $ and $ M\in\mathbb{R}^{(N_rN_\mu+1)\times N_rN_\mu} $ we use the pseudo-inverse yielding $ F = (M^tM)^{-1}M^tB $. Finally the use of a sparse matrix routine greatly improves the computation time. Our implementation was made using Matlab. The code is publicly available on GitLab [12].

    In our test cases we use the following parameters: $ r_{\min} = 0.06 $, $ r_{\max} = 5 $, $ \mu_{\min} = 0.05 $, $ \mu_{\max} = 5 $, $ \delta = 8 $, $ N_r = 70 $, $ N_\mu = 200 $, $ \kappa = 1 $, $ \nu = 1 $.

    On Figure 3 we compare the distribution functions $ f(r, \mu) $ obtained for various sets of parameters $ (p, \gamma) $. We also draw the corresponding marginal $ \rho(r) $ as well as $ \rho_0 $ and $ \rho_\text{fast} $. We observe that for small values of $ p $, $ \rho_\text{fast} $ is a good approximation of $ \rho $. For larger values it tends to amplify the phenomenon of variance reduction.

    Figure 3.  Numerical results. Numerical solution of the main Fokker-Planck model for various sets of parameters $ (\gamma, p) $. Left: Surface and contour plot of the distribution function $ f(r, \mu) $. The truncation at $ r = 2 $ and $ \mu = 2 $ is only for visualization purposes. Right: Corresponding marginal density $ \rho $ compared with $ \rho_\text{fast} $ and $ \rho_0 $.

    In order to confirm that the main Fokker-Planck model reduces the coefficient of variation as soon as $ \gamma > 0 $ we draw on Figure 4 the coefficient of variation for each distribution $ \rho, \rho_\text{fast} $ relatively to that of $ \rho_0 $ for several values of $ p $. We observe that indeed, the coefficient of variation is reduced. As in the case of fast $ \mu $RNA, the decay is more pronounced when the production of $ \mu $RNA is higher than that of mRNA, namely when $ p > 1 $. Interestingly enough, one also notices that the approximation $ \rho_\text{fast} $ increases the reduction of CV when $ p > 1 $ and diminishes it when $ p < 1 $. A transition at the special value $ p = 1 $ was already observed on Figure 1.

    Figure 4.  Numerical results. Relative coefficient of variation versus $ \gamma $ for various values of $ p $.

    In this section, we discuss the influence of the type of noise in the Fokker-Planck models. Let us go back to the system of stochastic differential equations considered at the beginning and generalize it as follows

    $ \left\{ drt=(crcrtμtkrrt)dt+ 2σrD(rt)dB1t,dμt=(cμcrtμtkμμt)dt+2σμD(μt)dB2t, \right. $

    with $ D $ some given function. In the models of the previous sections we chose $ D(x) = x^2 $. On the one hand it is natural to impose that $ D(x) $ vanishes when $ x\to0 $ in order to preserve the non-negativity of $ r_t $ and $ \mu_t $. On the other hand it is clear that the growth at infinity influences the tail of the equilibrium distribution which solves the corresponding Fokker-Planck equation. With a quadratic $ D $ we obtained algebraically decaying distributions.

    Nevertheless one may wonder if the decay of cell to cell variation due to $ \mu $RNA would still be observed if $ D $ is changed so that it involves distributions with faster decay at infinity. In order to answer this question, we choose a simple enough function $ D $ so that we can still derive analytical formulas for distributions of mRNA without binding and mRNA in the presence of "fast" $ \mu $RNA. Let us assume that

    $ D(r) = r\, . $

    In terms of modeling we may argue as in Section 2 and Section 3 in order to introduce the stationary probability distribution of mRNA without binding $ \tilde{\rho}_0 $ which solves

    $ \partial_r\left[\sigma_r\partial_r (r{\tilde{\rho}}_0) -(c_r - k_r\, r){\tilde{\rho}}_0\right]\ = \ 0\, . $

    It may still be solved analytically and one finds a gamma distribution

    $ ˜ρ0(r) = γcrσr,krσr(r) = Ccrσr,krσrrcrσr1ekrσrr $ (6.1)

    instead of an inverse gamma distribution in the quadratic case. The normalization constant is given in Section 3.1.

    In the case of fast $ \mu $RNA, we may once again follow the method of Section 2 and Section 3 and introduce $ {\tilde{\rho}}fast $ solving

    $ \partial_r\left[\partial_r(\sigma_r\ r\ {\tilde{\rho}}fast) - (\ c_r - c\ r\ {\tilde{j}}_\text{fast}(r) - k_r\ r){\tilde{\rho}}fast \right]\ = \ 0 $

    where the conditional expectation of the number of $ \mu $RNA within the population with $ r $ mRNA is given by

    $ {\tilde{j}}_\text{fast}(r)\ = \ \int_0^\infty \mu\, \gamma_{\frac{c_\mu}{\sigma_\mu}, \frac{k_\mu + c\, r}{\sigma_\mu}}(\mu)\, \mathrm{d}\mu\ = \ \frac{c_\mu}{k_\mu + c\, r}\, . $

    A direct computation then yields

    $ ˜ρfast(r) = C(1+ckμr)cμσrrcrσr1ekrσrr, $ (6.2)

    with $ C\equiv C(c_r, c_\mu, c, k_r, k\mu, \sigma_r, \sigma_\mu) $ a normalizing constant.

    Remark 6.1. Observe that the conditional expectation of the number of $ \mu $RNA within the population with $ r $ mRNA is unchanged, namely $ {\tilde{j}}_\text{fast}(r) = j_\text{fast}(r) $. More generally, the expectation of a univariate process $ (X_t)_t $ satisfying an SDE with linear drift $ \, \mathrm{d} X_t = (a + bX_t)\, \mathrm{d} t + \sqrt{2\sigma(X_t)}\, \mathrm{d} B_t $ does not depend on the diffusion coefficient $ \sigma $ as its density $ g $ satisfy $ \partial_t g(t, x) + \partial_x((a+bx)g(t, x)) - \partial^2_{xx}(\sigma(x)g(t, x)) = 0\, , $ so that multiplying by $ x $ and integrating yields $ \, \mathrm{d} \mathrm{E}[X_t] = (a + b\mathrm{E}[X_t])\, \mathrm{d} t $ on its expectation $ \mathrm{E}[X_t] $. The argument also holds for multivariate processes.

    Once again we seek the parameters of importance among the many parameters of the model by a dimensional analysis. The characteristic value of $ r $ remains $ \bar{r} = c_r/ k_r $ as it is the expectation of $ {\tilde{\rho}}_0 $. After rescaling we find the new distribution

    $ ˜ρη0(r) = γη,η(r) = Cη,ηrη1eηr, $ (6.3)

    and

    $ ˜ρfastη,γ,p(r) = Cadfastrη1(1+γr)pηeηr. $ (6.4)

    where the parameters $ p $ and $ \gamma $ are given by Eqs (4.7) and (4.6) respectively and still quantify the intensity of the binding and the respective production of $ \mu $RNA versus mRNA. The new parameter $ \eta $ is given by

    $ η = crσr. $ (6.5)

    In the context of a dimensional analysis, let us mention that it would be inaccurate to compare $ \eta $ and $ \delta $ as the $ \sigma_r $ (and $ \sigma_\mu $) do not represent the same quantity depending on the choice of $ D $. For $ D(r) = r^2 $ it has the same dimension as $ k_r $ so $ \delta = k_r/\sigma_r $ is the right dimensionless parameter. Here it has the same dimension as $ c_r $, which justifies the introduction of $ \eta $.

    The expectation, variance and coefficient of variation of $ {\tilde{\rho}}_0 $ are explicitly given by

    $ Exp(˜ρη0) = 1,Var(˜ρη0) = 1η,CV(˜ρη0) = 1η, $ (6.6)

    As there is no explicit formula for the coefficient of variation of $ {\tilde{\rho}}fast^{\eta, \gamma, p} $ we evaluate it numerically as in Section 4.3. The results are displayed on Figure 5. We observe that unlike the case of a quadratic diffusion coefficient the relative cell to cell variation, i.e., the cell to cell variation in the presence of $ \mu $RNA relative to cell to cell variations of the free case, is not unconditionally less than $ 1 $. For a large enough production of $ \mu $RNA, it eventually decays when the binding effect is very strong. However for smaller production of $ \mu $RNA or when the binding is weak, the effect is the opposite as the relative cell to cell variation is greater than $ 1 $. This is not satisfactory from the modeling point of view.

    Figure 5.  {Numerical computation of the cell to cell variation.} Relative cell to cell variation $ \mathrm{CV}({\tilde{\rho}}fast^{\eta, \gamma, p}) / \mathrm{CV}(\rho_0^{\eta}) $ for various parameters $ p $, $ \gamma $ and $ \eta $. On the horizontal axis, left means more production of mRNA and right means more production of $ \mu $RNA; On the vertical axis, top means more destruction of mRNA by binding and bottom means more destruction/consumption of mRNA by other mechanisms.

    In conclusion the choice of noise is important in this model. An unconditional cell to cell variation decay in the presence of $ \mu $RNA is observed for quadratic noise only. While other choices of noise may still lead to similar qualitative results, the choice $ D(r) = r^2 $ allowed us to derive explicit formulas for the approximate density $ \rho_\text{fast} $ which, as numerical simulations show, is fairly close to the marginal $ \rho $ corresponding to the solution of the main Fokker-Planck model.

    In this paper, we introduced a new model describing the joint probability density of the number of mRNA and $ \mu $RNA in a cell. It is based on a Fokker-Planck equation arising from a system of chemical kinetic equations for the number of two RNAs. The purpose of this simple model was to provide a mathematical framework to investigate how robustness in gene expression in a cell is affected by the presence of a regulatory feed-forward loop due to production of $ \mu $RNAs which bind to and deactivate mRNAs.

    Thanks to the combined use of analytical formulas and numerical simulations, we showed that robustness of gene expression is indeed affected by the presence of a feed-forward loop involving $ \mu $RNA production. However, whether the effect is regulatory or de-regulatory strongly depends on the assumptions made on the type of noise affecting both mRNA and $ \mu $RNAs production. In the case of geometric noise (the diffusivity being quadratic in the solution itself), the effect is to reduce the spread of the distribution as the reduction of the coefficient of variation shows. In the case of sub-geometric noise (the diffusivity being only linear in the solution itself), the effect increases the spread as shown by the increase of the coefficient of variation. We may attempt an explanation by comparing the mRNA distribution in the absence of $ \mu $RNA and in the limit of fast $ \mu $RNA in the two cases. In the quadratic diffusivity case, both distributions are fat-tailed (i.e., they decay polynomially with the number of unbound mRNA molecules $ r $, see Eqs (3.6) and (3.9)) but the rate of decay at infinity is modified by the presence of $ \mu $RNAs. On the other hand, in the linear diffusivity case, both decay exponentially fast (see Eqs (6.1) and (6.2)) and the exponential rate of decay is the same with or without $ \mu $RNAs. We propose that this might be the reason of the difference: In the quadratic diffusivity case, the change in polynomial decay allows to greatly reduce the standard deviation without affecting too much the mean, which results in a reduction of the coefficient of variation. In the linear diffusivity case, the exponential tail is not modified, which implies that the core of the distribution must be globally translated towards the origin, which affects the mean and the standard deviation in a similar way and does not systematically reduce the coefficient of variation. Which type of noise corresponds to the actual data is unknown at this stage. While quadratic diffusivity seems a fairly reasonable assumption (it is used in a number of contexts such as finance), it would require further experimental investigations to be fully justified in the present context. This discussion shows that the effect of $ \mu $RNA on noise regulation of mRNA translation is subtle and not easily predictable.

    Along the way we provided theoretical tools for the analysis of the Fokker Planck equation at play and robust numerical methods for simulations. As the main biological hypothesis for the usefulness of $ \mu $RNA in the regulation of gene expression is based on their ability to reduce external noise, we also discussed the particular choice of stochasticity in the model.

    There are several perspectives to this work. A first one would be the calibration of the parameters of the model from real-world data. This would allow to quantify more precisely the amount of cell-to-cell variation reduction due to $ \mu $RNA, thanks to the thorough numerical investigation done in this contribution of the effects of the parameters of the model. Besides, another perspective would be an improvement of inequality (4.12) to the conjecture (4.3). This would bring a definitive theoretical answer to the hypothesis of increased gene expression level in the simplified model of "fast" $ \mu $RNAs. One may also look into establishing a similar inequality for the general model. Finally, the gene regulatory network in a cell is considerably more complex than the simple, yet enlightening in our opinion, dynamics proposed in this paper. A natural improvement would be the consideration of more effects in the model, such as the production of the transcription factor, or the translation of mRNA into proteins, among many others.

    PD acknowledges support by the Engineering and Physical Sciences Research Council (EPSRC) under grants no. EP/M006883/1 and EP/N014529/1, by the Royal Society and the Wolfson Foundation through a Royal Society Wolfson Research Merit Award no. WM130048 and by the National Science Foundation (NSF) under grant no. RNMS11-07444 (KI-Net). PD is on leave from CNRS, Institut de Mathématiques de Toulouse, France. MH acknowledges support by the Labex CEMPI (ANR-11-LABX-0007-01). SM acknowledges support by the CNRS–Royal Society exchange projects "CODYN" and "Segregation models in social sciences" and the Chaire Modélisation Mathématique et Biodiversité of Véolia Environment - École Polytechnique - Museum National d’Histoire Naturelle - Fondation X. All authors would like to thank Prof. Matthias Merkenshlager from Imperial College Institute of Clinical Sciences for bringing this problem to their attention and stimulating discussions.

    The authors declare no conflict of interest.

    No new data were collected in the course of this research.

    In this section we give a elementary proof of the 1D version of the Brascamp and Lieb inequality (see [23,Theorem 4.1], [24]), which is an extension of the Gaussian Poincaré inequality in the case of log-concave measures. This allows us to derive a weighted Poincaré inequality for the gamma distribution and deduce, by a change of variable, a similar functional inequality for the inverse gamma distribution.

    Proposition A.1. Let $ I\subset\mathbb{R} $ be an open non-empty interval and $ V:I\rightarrow\mathbb{R} $ a function of class $ \mathcal{C}^2 $. Assume that

    (i) $ V $ is strictly convex;

    (ii) $ e^{-V} $ is a probability density on $ I $;

    (iii) $ V $ tends to $ +\infty $ at the extremities of $ I $.

    Then, for any suitably integrable function $ u $, one has

    $ I|u(x)ueV|2eV(x)dx  I|u(x)|2eV(x)(V"(x))1dx, $ (A.1)

    where for a density $ \nu $ the notation $ \langle\, u\, \rangle_{\nu} $ denotes $ \int u\nu $.

    Proof. Without loss of generality, as one may replace $ u $ with $ u-\langle\, u\, \rangle_{e^{-V}} $, we assume that $ \langle\, u\, \rangle_{e^{-V}} = 0 $. We also assume that $ u $ is of class $ \mathcal{C}^1 $ and compactly supported in $ I $ and one can then extend a posteriori the class of admissible function by a standard density argument. Then, using (ii) one has

    $ I|u(x)|2eV(x)dx=12I×I|u(x)u(y)|2e(V(x)+V(y))dxdy=12I×I|yxu(z)dz|2e(V(x)+V(y))dxdy. $

    Now using the Cauchy-Schwarz inequality and assumption (i) one has

    $ \int_I |u(x)|^2\, e^{-V(x)}\, \, \mathrm{d} x\ \leq \frac{1}{2}\, \iint_{I\times I}\left(\int_x^y|u'(z)|^2 (V"(z))^{-1}\, \, \mathrm{d} z\right)\left(\int_x^yV"(z)\, \mathrm{d} z\right)\, e^{-(V(x)+V(y))}\, \, \mathrm{d} x\, \, \mathrm{d} y\, . $

    Then take any point $ x_0\in I $ and define

    $ U(x) = \int_{x_0}^x|u'(z)|^2\, (V"(z))^{-1}\, \, \mathrm{d} z\, , $

    so that the inequality rewrites

    $ \int_I |u(x)|^2\, e^{-V(x)}\, \, \mathrm{d} x \leq \frac{1}{2}\, \iint_{I\times I}\left(U(y)-U(x)\right)\left(V'(y)-V'(x)\right)\, e^{-(V(x)+V(y))}\, \, \mathrm{d} x\, \, \mathrm{d} y\, . $

    Now just expand the right-hand side and use Fubini's theorem on each term as well as assumptions (ii) and (iii) to obtain

    $ \int_I |u(x)|^2\, e^{-V(x)}\, \, \mathrm{d} x\ \leq\ \int_I U(x)\, V'(x)\, e^{-V(x)}\, \, \mathrm{d} x\, . $

    One concludes by integrating the right-hand side by parts and observing that boundary terms vanish again by assumption (iii).

    Remark A.2. The proof is an adaptation of the original proof of the (flat) Poincaré-Wirtinger inequality by Poincaré.

    Observe that for $ I = \mathbb{R} $ and $ V(x) = x^2/2 $, one recovers the classical Gaussian Poincaré inequality.

    Remark A.3. The inequality is sharp. It is an equality for functions of the form $ u(x) = aV'(x)+b $, with $ a, b\in\mathbb{R} $ if $ V $ is such that $ V'(x)e^{-V(x)} $ tends to $ 0 $ at the boundaries, and only for constant functions otherwise (i.e., $ a = 0 $ and $ b\in\mathbb{R} $).

    From the Brascamp-Lieb inequality, we now infer Poincaré inequalities for gamma and inverse gamma distributions.

    Proposition A.4. Let $ \alpha > 1 $ and $ \beta > 0 $. Then, for any functions $ u, v $ such that the integrals make sense, one has

    $ 0|u(x)uγα,β|2γα,β(x)dx  1α10|u(x)|2γα,β(x)x2dx, $ (A.2)

    and

    $ 0|v(y)vgα,β|2gα,β(y)dy  1α10|v(y)|2gα,β(y)y2dy, $ (A.3)

    where for a probability density $ \nu $ the notation $ \langle\, u\, \rangle_{\nu} $ denotes $ \int u\nu $.

    Proof. The first inequality is an application of Eq (A.1) with $ V(x) = \beta x - (\alpha-1)\ln(x) - \ln(C_{\alpha, \beta}) $, where $ C_{\alpha, \beta} = \beta^\alpha/\Gamma(\alpha) $. Then take $ v(y) = u(1/y) $ and make the change of variable $ y = 1/x $ in all the integrals of Eq (A.2) to get the result.

    Remark A.5. To the best of our knowledge the classical Bakry and Emery method does not seem to apply to show directly the functional inequalities of Proposition A.4. Let us give some details. In order to show a Poincaré inequality of the type

    $ \int_I|u(x) - \langle\, u\, \rangle_{e^{-V}}|^2\, e^{-V(x)}\, \, \mathrm{d} x\ \leq\ \int_I|u'(x)|^2\, e^{-V(x)}\, D(x)\, \, \mathrm{d} x\, , $

    for $ V $ as in Proposition A.1, it is sufficient that $ D $ and $ V $ satisfy the following curvature-dimension inequality

    $ R(x) := 14(D(x))212D"(x)D(x)+D(x)2V"(x)+12D(x)D(x)V(x)  λ1D(x), $ (A.4)

    for some positive constant $ \lambda_1 > 0 $. We refer to [25,26] for the general form of the latter Bakry-Emery condition (for multidimensional anisotropic inhomogeneous diffusions) and to [27] or [28] for the simpler expression in the case of isotropic inhomogeneous diffusion, as discussed here. In the case of the inequalities (A.2) and (A.3), one has respectively $ D(x) = x^2/(\alpha-1) $, $ V(x) = \beta x - (\alpha-1)\ln(x) - \ln(C_{\alpha, \beta}) $ and $ D(y) = y^2/(\alpha-1) $, $ V(y) = \beta / y + (\alpha+1)\ln(y) - \ln(C_{\alpha, \beta}) $, which yields respectively $ R(x) = \beta\, x^3/(\alpha-1)^2 $ and $ R(y) = \beta\, y/(\alpha-1)^2 $. As claimed above, neither Eq (A.2) nor Eq (A.3) satisfy the condition (A.4). One also observes that in both cases the curvature-dimension inequality fails because of a degeneracy at one end of the interval.

    Let us finally mention that there are in the literature other occurrences of Poincaré and more generally convex Sobolev inequalities for gamma distributions [29,30,31,32]. However, we found out that the diffusion coefficient is always taken of the form $ D(x) = 4x/\beta $. This weight, associated with the gamma invariant measure, corresponds to the Laguerre diffusion $ L_{\alpha, \beta}f (x) = \beta x\, f" (\beta x) - (\alpha - \beta x)f'(\beta x) $. This operator differs from the adjoint of the one appearing in our model given by Eq (2.6). In this case, one can check that the curvature-dimension condition of Bakry and Emery is satisfied as soon as $ \alpha\geq 1/2 $.

    Proposition A.6. One has the bound

    $ \mathrm{CV}(\rho_\mathit{\text{fast}}^{\delta, \gamma, p})\ \leq\ C_\delta\ : = \ \left(\left(\frac{\delta}{\delta-1}\right)^2\left(1-\frac{1}{(\delta-1)^2}\right)^{{\delta-2}}-1\right)^{\frac12}\, , $

    which holds for all $ \delta > 2 $, $ \gamma > 0 $ and $ p\geq0 $.

    Proof. The bound is a consequence of the Prékopa-Leindler inequality (see [23] and references therein) which states that if $ f, g, h:\mathbb{R}^d\to[0, +\infty) $ are three functions satisfying for some $ \lambda\in(0, 1) $ and for all $ x, y $,

    $ h((1λ)x+λy)f(x)1λg(y)λ, $ (A.5)

    then

    $ hL1(Rd)  f1λL1(Rd)gλL1(Rd). $ (A.6)

    We use it with $ \lambda = 1/2 $, $ f(x) = (1+x)^{\gamma p\delta}x^{\delta-2}e^{-\delta \gamma x} $ if $ x\geq0 $ and $ f(x) = 0 $ if $ x < 0 $, $ g(x) = x^2f(x) $ and $ h(x) = (1+C_\delta^2)^{1/2}xf(x) $. The condition (A.5) is then equivalent to

    $ (1+C_\delta^2)^{-1/2}\ \leq\ \left[\frac{\left(1+\frac{x+y}{2}\right)}{(1+x)(1+y)}\right]^{\frac{\gamma p\delta}{2}}\left(\frac yx\right)^{\frac12}\left(\frac{\left(\frac yx\right)^{\frac12}+\left(\frac xy\right)^{\frac12}}{2}\right)^{\delta-1} $

    which is satisfied as the term between brackets is always greater than $ 1 $ and the function $ z\mapsto z[(z+z^{-1})/2]^{\delta-1} $, $ z > 0 $ is bounded from below by $ (1+C_\delta^2)^{-1/2} $, where $ C_\delta $ is given in Eq (4.12). Then with the change of variable $ x' = 1/(\gamma x) $ in the integrals of Eq (A.6), one recovers Eq (4.12).

    [1] James W (1890) The principles of psychology (H. Holt, New York, NY).
    [2] Posner MI, Petersen SE (1990) The attention system of the human brain. Annu Rev Neurosci 13: 25-42. doi: 10.1146/annurev.ne.13.030190.000325
    [3] Norman DA, Shallice T (1986) Attention to action: willed and automatic control of behavior. Consciousness and Self-Regulation, eds Davison RJ, Schwartz GE, Shapiro D (Plenum Press, New York, NY), pp 1-18.
    [4] Corbetta M, Shulman GL (2002) Control of goal-directed and stimulus-driven attention in the brain. Nat Rev Neurosci 3: 201-215.
    [5] Posner MI, DiGirolamo GJ (1998) Executive attention: Conflict, target detection, and cognitive control. The Attentive Brain, ed Parasuraman R (MIT Press, Cambridge, MA), pp 401-423.
    [6] D’Angelo MC, Milliken B, Jiménez L, et al. (2013) Implementing flexibility in automaticity: Evidence from context-specific implicit sequence learning. Conscious Cogn 22: 64-81. doi: 10.1016/j.concog.2012.11.002
    [7] Rueda MR, Posner MI, Rothbart MK (2005) The development of executive attention: contributions to the emergence of self-regulation. Dev Neuropsychol 28: 573-94. doi: 10.1207/s15326942dn2802_2
    [8] Luu P, Tucker DM, Derryberry D, et al. (2003) Electrophysiological responses to errors and feedback in the process of action regulation. Psychol Sci 14: 47-53. doi: 10.1111/1467-9280.01417
    [9] Dagenbach D, Carr TH (1994) Inhibitory processes in attention, memory, and language (Academic Press, San Diego, CA).
    [10] Petersen SE, Posner MI (2012) The Attention System of the Human Brain: 20 Years After. Annu Rev Neurosci 35: 73-89. doi: 10.1146/annurev-neuro-062111-150525
    [11] Posner MI, Rueda MR, Kanske P (2007) Probing the Mechanisms of Attention. Handb Psychophysiol: 410-432.
    [12] Fan J, McCandliss BD, Sommer T, et al. (2002) Testing the efficiency and independence of attentional networks. J Cogn Neurosci 14: 340-347. doi: 10.1162/089892902317361886
    [13] Callejas A, Lupiàñez J, Funes MJ, et al. (2005) Modulations among the alerting, orienting and executive control networks. Exp brain Res 167: 27-37. doi: 10.1007/s00221-005-2365-z
    [14] Fan J, Gu X, Guise KG, et al. (2009) Testing the behavioral interaction and integration of attentional networks. Brain Cogn 70: 209-220. doi: 10.1016/j.bandc.2009.02.002
    [15] Hackley SA, Valle-Inclán F (2003) Which stages of processing are speeded by a warning signal? Biol Psychol 64: 27-45. doi: 10.1016/S0301-0511(03)00101-7
    [16] Weinbach N, Henik A (2012) The relationship between alertness and executive control. J Exp Psychol Hum Percept Perform 38: 1530-1540. doi: 10.1037/a0027875
    [17] Pozuelos JP, Paz-Alonso PM, Castillo A, et al. (2014) Development of Attention Networks and Their Interactions in Childhood. Dev Psychol 50: 2405-2415. doi: 10.1037/a0037469
    [18] Fox MD, Corbetta M, Snyder AZ, et al. (2006) Spontaneous neuronal activity distinguishes human dorsal and ventral attention systems. Proc Natl Acad Sci U S A 103: 10046-10051. doi: 10.1073/pnas.0604187103
    [19] Dosenbach NUF, Fair DA, Miezin FM, et al. (2007) Distinct brain networks for adaptive and stable task control in humans. Proc Natl Acad Sci U S A 104: 11073-11078. doi: 10.1073/pnas.0704320104
    [20] Posner MI, Rothbart MK (2007) Research on attention networks as a model for the integration of psychological science. Annu Rev Psychol 58: 1-23. doi: 10.1146/annurev.psych.58.110405.085516
    [21] Coull JT, Nobre AC, Frith CD (2001) The noradrenergic a2 agonist clonidine modulates behavioural and neuroanatomical correlates of human attentional orienting and alerting. Cereb Cortex 11: 73-84. doi: 10.1093/cercor/11.1.73
    [22] Aston-Jones G, Cohen JD (2005) An integrative theory of locus coeruleus-norepinephrine function: adaptive gain and optimal performance. Annu Rev Neurosci 28: 403-450. doi: 10.1146/annurev.neuro.28.061604.135709
    [23] Coull JT, Frith CD, Frackowiak RSJ, Grasby PM (1996) A fronto-parietal network for rapid visual information processing: A PET study of sustained attention and working memory. Neuropsychologia 34: 1085-1095. doi: 10.1016/0028-3932(96)00029-2
    [24] Cui RQ, Egkher A, Huter D, et al. (2000) High resolution spatiotemporal analysis of the contingent negative variation in simple or complex motor tasks and a non-motor task. Clin Neurophysiol 111: 1847-1859. doi: 10.1016/S1388-2457(00)00388-6
    [25] Coull JT (2004) fMRI studies of temporal attention: Allocating attention within, or towards, time. Cogn Brain Res 21: 216-226. doi: 10.1016/j.cogbrainres.2004.02.011
    [26] Hillyard SA (1985) Electrophysiology of human selective attention. Trends Neurosci 8: 400-405. doi: 10.1016/0166-2236(85)90142-0
    [27] Mangun GR, Hillyard SA (1987) The spatial allocation of visual attention as indexed by event-related brain potentials. Hum Factors 29: 195-211.
    [28] Desimone R, Duncan J (1995) Neural mechanisms of selective visual attention. Annu Rev Neurosci 18: 193-222. doi: 10.1146/annurev.ne.18.030195.001205
    [29] Corbetta M, Patel G, Shulman GL (2008) The Reorienting System of the Human Brain: From Environment to Theory of Mind. Neuron 58 : 306-324.
    [30] Greicius MD, Krasnow B, Reiss AL, et al. (2003) Functional connectivity in the resting brain: a network analysis of the default mode hypothesis. Proc Natl Acad Sci U S A 100: 253-258. doi: 10.1073/pnas.0135058100
    [31] Mantini D, Perrucci MG, Del Gratta C, et al. (2007) Electrophysiological signatures of resting state networks in the human brain. Proc Natl Acad Sci U S A 104: 13170-13175. doi: 10.1073/pnas.0700668104
    [32] Visintin E, De Panfilis C, Antonucci C, et al. (2015) Parsing the intrinsic networks underlying attention: A resting state study. Behav Brain Res 278: 315-322. doi: 10.1016/j.bbr.2014.10.002
    [33] Umarova RM, Saur D, Schnell S, et al. (2010) Structural connectivity for visuospatial attention: Significance of ventral pathways. Cereb Cortex 20: 121-129. doi: 10.1093/cercor/bhp086
    [34] Buschman TJ, Miller EK (2007) Top-Down Versus Bottom-Up Control of Attention in the Prefrontal and Posterior Parietal Cortices. Sci 315: 1860-1862. doi: 10.1126/science.1138071
    [35] Vossel S, Geng JJ, Fink GR (2013) Dorsal and Ventral Attention Systems: Distinct Neural Circuits but Collaborative Roles. Neurosci 20: 150-159.
    [36] He BJ, AZ Snyder, JL Vincent, et al. (2007) Breakdown of functional connectivity in frontoparietal networks underlies behavioral deficits in spatial neglect. Neuron 53: 905-918. doi: 10.1016/j.neuron.2007.02.013
    [37] Giesbrecht B, Weissman DH, Woldorff MG, et al. (2006) Pre-target activity in visual cortex predicts behavioral performance on spatial and feature attention tasks. Brain Res 1080: 63-72. doi: 10.1016/j.brainres.2005.09.068
    [38] Geng JJ, Mangun GR (2011) Right temporoparietal junction activation by a salient contextual cue facilitates target discrimination. Neuroimage 54: 594-601. doi: 10.1016/j.neuroimage.2010.08.025
    [39] Fan J, Flombaum JI, McCandliss BD, et al. (2003) Cognitive and Brain Consequences of Conflict. Neuroimage 18: 42-57. doi: 10.1006/nimg.2002.1319
    [40] Bush G, Luu P, Posner MI (2000) Cognitive and emotional influences in anterior cingulate cortex. Trends Cogn Sci 4: 215-222. doi: 10.1016/S1364-6613(00)01483-2
    [41] Drevets WC, Raichle ME (1998) Reciprocal suppression of regional cerebral blood flow during emotional versus higher cognitive processes: Implications for interactions between emotion and cognition. Cogn emotioin 12: 353-385. doi: 10.1080/026999398379646
    [42] Botvinick MM, Nystrom L, Fissell K, et al. (1999) Conflict monitoring versus selection-for-action in anterior cingulate cortex. Nature 402: 179-181. doi: 10.1038/46035
    [43] Botvinick MM, Braver TS, Barch DM, et al. (2001) Conflict monitoring and cognitive control. Psychol Rev 108: 624-652. doi: 10.1037/0033-295X.108.3.624
    [44] Kopp B, Tabeling S, Moschner C, et al. (2006) Fractionating the Neural Mechanisms of Cognitive Control. J Cogn Neurosci: 949-965.
    [45] Van Veen V, Carter CS (2002) The timing of action-monitoring processes in the anterior cingulate cortex. J Cogn Neurosci 14: 593-602. doi: 10.1162/08989290260045837
    [46] Posner MI, Sheese BE, Odludaş Y, et al. (2006) Analyzing and shaping human attentional networks. Neural Networks 19: 1422-1429. doi: 10.1016/j.neunet.2006.08.004
    [47] Dosenbach NUF, Fair Da, Cohen AL, et al. (2008) A dual-networks architecture of top-down control. Trends Cogn Sci 12: 99-105. doi: 10.1016/j.tics.2008.01.001
    [48] Rueda MR (2014) Development of Attention. Oxford Handb Cogn Neurosci 1: 296-318.
    [49] Rueda MR, Posner MI, Rothbart MK, et al. (2004) Development of the time course for processing conflict: an event-related potentials study with 4 year olds and adults. BMC Neurosci 5: 39. doi: 10.1186/1471-2202-5-39
    [50] Abundis-Gutiérrez A, Checa P, Castellanos C, et al. (2014) Electrophysiological correlates of attention networks in childhood and early adulthood. Neuropsychologia 57: 78-92. doi: 10.1016/j.neuropsychologia.2014.02.013
    [51] Gießing C, Thiel CM, Alexander-Bloch aF, et al. (2013) Human brain functional network changes associated with enhanced and impaired attentional task performance. J Neurosci 33: 5903-5914. doi: 10.1523/JNEUROSCI.4854-12.2013
    [52] Gao W, Zhub HT, Giovanello KS, et al. (2009) Evidence on the emergence of the brain’s default network from 2-week-old to 2-year-old healthy pediatric subjects. Proc Natl Acad Sci U S A 106: 6790-6795. doi: 10.1073/pnas.0811221106
    [53] Fair DA, Dosenbach NU, Church JA, et al. (2007) Development of distinct control networks through segregation and integration. Proc Natl Acad Sci U S A 104: 13507-13512. doi: 10.1073/pnas.0705843104
    [54] Fair DA, Cohen AL, Power JD, et al. (2009) Functional brain networks develop from a “local to distributed” organization. PLoS Comput Biol 5: e1000381. doi: 10.1371/journal.pcbi.1000381
    [55] Fan J, Wu Y, Fossella JA, et al. (2001) Assessing the heritability of attentional networks. BMC Neurosci 2: 14. doi: 10.1186/1471-2202-2-14
    [56] Marrocco RT, Davidson MC (1998) Neurochemistry of attention. The Attentive Brain, ed Parasuraman R (MIT Press, Cambridge, MA), pp 35-50.
    [57] Congdon E, Lesch KP, Canli T (2008) Analysis of DRD4 and DAT polymorphisms and behavioral inhibition in healthy adults: implications for impulsivity. Am J Med Genet B Neuropsychiatr Genet 147: 27-32.
    [58] Rueda MR, Rothbart MK, McCandliss BD, et al. (2005) Training, maturation, and genetic influences on the development of executive attention. Proc Natl Acad Sci U S A 102: 14931-14936. doi: 10.1073/pnas.0506897102
    [59] Diamond A (2007) Consequences of variations in genes that affect dopamine in prefrontal cortex. Cereb cortex 17: i161-170. doi: 10.1093/cercor/bhm082
    [60] Forbes EE, Brown SM, Kimak M, et al. (2009) Genetic variation in components of dopamine neurotransmission impacts ventral striatal reactivity associated with impulsivity. Mol Psychiatry 14: 60-70. doi: 10.1038/sj.mp.4002086
    [61] Congdon E, Constable RT, Lesch KP, et al. (2009) Influence of SLC6A3 and COMT variation on neural activation during response inhibition. Biol Psychol 81: 144-152. doi: 10.1016/j.biopsycho.2009.03.005
    [62] Mueller EM, Makeig S, Stemmler G, et al. (2011) Dopamine effects on human error processing depend on catechol-O-methyltransferase VAL158MET genotype. J Neurosci 31: 15818-15825. doi: 10.1523/JNEUROSCI.2103-11.2011
    [63] Espeseth T, Sneve MH, Rootwelt H, et al. (2010) Nicotinic receptor gene CHRNA4 interacts with processing load in attention. PLoS One 5: e14407. doi: 10.1371/journal.pone.0014407
    [64] Greenwood PM, Parasuraman R, Espeseth T (2012) A cognitive phenotype for a polymorphism in the nicotinic receptor gene CHRNA4. Neurosci Biobeha Rev 36: 1331-1341. doi: 10.1016/j.neubiorev.2012.02.010
    [65] Lundwall Ra, Guo DC, Dannemiller JL (2012) Exogenous visual orienting is associated with specific neurotransmitter genetic markers: A population-based genetic association study. PLoS One 7.
    [66] Zozulinsky P, Greenbaum L, Brande-Eilat N, et al. (2014) Dopamine system genes are associated with orienting bias among healthy individuals. Neuropsychologia 62: 48-54. doi: 10.1016/j.neuropsychologia.2014.07.005
    [67] Sheese BE, Voelker P, Posner MI, et al. (2009) Genetic variation influences on the early development of reactive emotions and their regulation by attention. Cogn Neuropsychiatry 14: 332-355. doi: 10.1080/13546800902844064
    [68] Posner MI, Rothbart MK, Sheese BE (2007) Attention genes. Dev Sci 10: 24-29. doi: 10.1111/j.1467-7687.2007.00559.x
    [69] Bornstein MH, Bradley RH (2003) Socioeconomic status, parenting, and child development (Lawrence Erlbaum Associates Publishers, Mahwah, NJ).
    [70] Bernier A, Carlson SM, Whipple N (2010) From external regulation to self-regulation: early parenting precursors of young children’s executive functioning. Child Dev 81: 326-339. doi: 10.1111/j.1467-8624.2009.01397.x
    [71] Gaertner BM, Spinrad TL, Eisenberg N (2008) Focused attention in toddlers: Measurement, stability, and relations to negative emotion and parenting. Infant Child Dev 17: 339-363. doi: 10.1002/icd.580
    [72] Cipriano EA, Stifter CA (2010) Predicting preschool effortful control from toddler temperament and parenting behavior. J Appl Dev Psychol 31: 221-230. doi: 10.1016/j.appdev.2010.02.004
    [73] Liew J, Chen Q, Hughes JN (2010) Child Effortful Control, Teacher-student Relationships, and Achievement in Academically At-risk Children: Additive and Interactive Effects. Early Child Res Q 25: 51-64. doi: 10.1016/j.ecresq.2009.07.005
    [74] Hackman D, Farah M (2009) Socioeconomic status and the developing brain Daniel. Trends Cogn Sci 13: 65-73. doi: 10.1016/j.tics.2008.11.003
    [75] Wanless SB, McClelland MM, Tominey SL, et al. (2011) The Influence of Demographic Risk Factors on Children’s Behavioral Regulation in Prekindergarten and Kindergarten. Early Educ Dev 22: 461-488. doi: 10.1080/10409289.2011.536132
    [76] Mezzacappa E (2004) Alerting, orienting, and executive attention: developmental properties and sociodemographic correlates in an epidemiological sample of young, urban children. Child Dev 75: 1373-1386. doi: 10.1111/j.1467-8624.2004.00746.x
    [77] Clearfield MW, Niman LC (2012) SES affects infant cognitive flexibility. Infant Behav Dev 35: 29-35. doi: 10.1016/j.infbeh.2011.09.007
    [78] Lawson GM, Duda JT, Avants BB, et al. (2013) Associations between children’s socioeconomic status and prefrontal cortical thickness. Dev Sci 16: 641-652. doi: 10.1111/desc.12096
    [79] Jolles DD, Crone EA (2012) Training the developing brain: a neurocognitive perspective. Front Hum Neurosci 6: 76.
    [80] Tang Y-Y, Posner MI (2009) Attention training and attention state training. Trends Cogn Sci 13: 222-227. doi: 10.1016/j.tics.2009.01.009
    [81] Karbach J, Kray J (2009) How useful is executive control training? Age differences in near and far transfer of task-switching training. Dev Sci 12: 978-990.
    [82] Jaeggi SM, Buschkuehl M, Jonides J, et al. (2011) Short- and long-term benefits of cognitive training. Proc Natl Acad Sci U S A 108: 10081-10086. doi: 10.1073/pnas.1103228108
    [83] Thorell LB, Lindqvist S, Bergman Nutley S, et al. (2008) Training and transfer effects of executive functions in preschool children. Dev Sci 12: 106-113.
    [84] Olesen PJ, Westerberg H, Klingberg T (2004) Increased prefrontal and parietal activity after training of working memory. Nat Neurosci 7: 75-79. doi: 10.1038/nn1165
    [85] Jolles DD, Van Buchem MA, Crone EA, et al. (2013) Functional brain connectivity at rest changes after working memory training. Hum Brain Mapp 34: 396-406. doi: 10.1002/hbm.21444
    [86] McNab F, Andrea V, Lars F, et al. (2009) Changes in cortical dopamine D1 receptor binding associated with cognitive training. Science 323: 800-802. doi: 10.1126/science.1166102
    [87] Tang Y-Y, Posner MI (2014) Training brain networks and states. Trends Cogn Sci 18: 345-350. doi: 10.1016/j.tics.2014.04.002
    [88] Malinowski P (2013) Neural mechanisms of attentional control in mindfulness meditation. Front Neurosci 7: 8.
    [89] Tang Y-Y, Ma YH, Wang JH, et al. (2007) Short-term meditation training improves attention and self-regulation. Proc Natl Acad Sci U S A 104: 17152-17156. doi: 10.1073/pnas.0707678104
    [90] Moore A, Gruber T, Derose J, et al. (2012) Regular, brief mindfulness meditation practice improves electrophysiological markers of attentional control. Front Hum Neurosci 6: 1-15.
    [91] Slagter HA, Lutz A, Greischar LL, et al. (2007) Mental Training Affects Distribution of Limited Brain Resources. Plos Biol 5.
    [92] Hölzel BK, Ott U, Hempel H, et al. (2007) Differential engagement of anterior cingulate and adjacent medial frontal cortex in adept meditators and non-meditators. Neurosci Lett 421: 16-21. doi: 10.1016/j.neulet.2007.04.074
    [93] Tang Y-Y, Lu Q, Fan M, et al. (2012) Mechanisms of white matter changes induced by meditation. Proc Natl Acad Sci: 1-5.
    [94] Posner MI, Tang Y-Y, Lynch G (2014) Mechanisms of white matter change induced by meditation training. Front Psychol 5: 1-4.
    [95] Tang Y-Y, Lua Ql, Gengc XJ, et al. (2010) Short-term meditation induces white matter changes in the anterior cingulate. Proc Natl Acad Sci U S A 107: 15649-15652. doi: 10.1073/pnas.1011043107
    [96] Hebb DO (1949) Organization of behavior (John Wiley & Sons, New York, NY).
  • This article has been cited by:

    1. Yekaterina Epshteyn, Chang Liu, Chun Liu, Masashi Mizuno, Nonlinear inhomogeneous Fokker–Planck models: Energetic-variational structures and long-time behavior, 2022, 20, 0219-5305, 1295, 10.1142/S0219530522400036
    2. Yekaterina Epshteyn, Chang Liu, Chun Liu, Masashi Mizuno, Local well-posedness of a nonlinear Fokker–Planck model, 2023, 36, 0951-7715, 1890, 10.1088/1361-6544/acb7c2
    3. Nicolas Boulanger, Fabien Buisseret, Victor Dehouck, Frédéric Dierick, Olivier White, Diffusion in Phase Space as a Tool to Assess Variability of Vertical Centre-of-Mass Motion during Long-Range Walking, 2023, 5, 2624-8174, 168, 10.3390/physics5010013
  • Reader Comments
  • © 2015 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(17373) PDF downloads(3047) Cited by(47)

Figures and Tables

Figures(5)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog