Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Stochastic simulations of the Schnakenberg model with spatial inhomogeneities using reactive multiparticle collision dynamics

  • A numerically efficient globally averaged number density approach is used to simulate a reaction-diffusion system using a particle-based stochastic simulation algorithm called reactive multiparticle collision (RMPC) dynamics. Constant diffusivity of the particles is achieved through a time-varying rotation angle (also called collision angle). Variation in the diffusion coefficient between two different chemical species is achieved in one of two ways: (ⅰ) using a different kBT/m value for one species compared to the other, or (ⅱ) using the same kBT/m value for both species, but using a different probability to free-stream for one species compared to another. For smaller diffusivities and larger spatial inhomogeneities, bath particles were necessary for the model to agree with the PDE solution. The latter approach was further used without a bath, and shown to be capable of producing Turing patterns after long simulation times. The significance of our work is that RMPC can serve as a feasible simulation tool for both short and long-term simulations, can handle spatial inhomogeneities, can model a fairly large range of diffusivities in a reaction-diffusion scenario, and is capable of producing Turing patterns. An advantage of this method includes more detailed system information in feasible simulation times.

    Citation: Alireza Sayyidmousavi, Katrin Rohlf. Stochastic simulations of the Schnakenberg model with spatial inhomogeneities using reactive multiparticle collision dynamics[J]. AIMS Mathematics, 2019, 4(6): 1805-1823. doi: 10.3934/math.2019.6.1805

    Related Papers:

    [1] Mhamed Eddahbi, Mogtaba Mohammed, Hammou El-Otmany . Well-posedness of a nonlinear stochastic model for a chemical reaction in porous media and applications. AIMS Mathematics, 2024, 9(9): 24860-24886. doi: 10.3934/math.20241211
    [2] Anna Sun, Ranchao Wu, Mengxin Chen . Turing-Hopf bifurcation analysis in a diffusive Gierer-Meinhardt model. AIMS Mathematics, 2021, 6(2): 1920-1942. doi: 10.3934/math.2021117
    [3] Yongchang Wei, Zongbin Yin . Long-time dynamics of a stochastic multimolecule oscillatory reaction model with Poisson jumps. AIMS Mathematics, 2022, 7(2): 2956-2972. doi: 10.3934/math.2022163
    [4] Xiaoming Wang, Muhammad W. Yasin, Nauman Ahmed, Muhammad Rafiq, Muhammad Abbas . Numerical approximations of stochastic Gray-Scott model with two novel schemes. AIMS Mathematics, 2023, 8(3): 5124-5147. doi: 10.3934/math.2023257
    [5] Jinji Du, Chuangliang Qin, Yuanxian Hui . Optimal control and analysis of a stochastic SEIR epidemic model with nonlinear incidence and treatment. AIMS Mathematics, 2024, 9(12): 33532-33550. doi: 10.3934/math.20241600
    [6] Tahir Khan, Fathalla A. Rihan, Muhammad Bilal Riaz, Mohamed Altanji, Abdullah A. Zaagan, Hijaz Ahmad . Stochastic epidemic model for the dynamics of novel coronavirus transmission. AIMS Mathematics, 2024, 9(5): 12433-12457. doi: 10.3934/math.2024608
    [7] Ruonan Liu, Tomás Caraballo . Random dynamics for a stochastic nonlocal reaction-diffusion equation with an energy functional. AIMS Mathematics, 2024, 9(4): 8020-8042. doi: 10.3934/math.2024390
    [8] Dan-Dan Dai, Wei Zhang, Yu-Lan Wang . Numerical simulation of the space fractional (3+1)-dimensional Gray-Scott models with the Riesz fractional derivative. AIMS Mathematics, 2022, 7(6): 10234-10244. doi: 10.3934/math.2022569
    [9] Heping Jiang . Complex dynamics induced by harvesting rate and delay in a diffusive Leslie-Gower predator-prey model. AIMS Mathematics, 2023, 8(9): 20718-20730. doi: 10.3934/math.20231056
    [10] Jawdat Alebraheem . Asymptotic stability of deterministic and stochastic prey-predator models with prey herd immigration. AIMS Mathematics, 2025, 10(3): 4620-4640. doi: 10.3934/math.2025214
  • A numerically efficient globally averaged number density approach is used to simulate a reaction-diffusion system using a particle-based stochastic simulation algorithm called reactive multiparticle collision (RMPC) dynamics. Constant diffusivity of the particles is achieved through a time-varying rotation angle (also called collision angle). Variation in the diffusion coefficient between two different chemical species is achieved in one of two ways: (ⅰ) using a different kBT/m value for one species compared to the other, or (ⅱ) using the same kBT/m value for both species, but using a different probability to free-stream for one species compared to another. For smaller diffusivities and larger spatial inhomogeneities, bath particles were necessary for the model to agree with the PDE solution. The latter approach was further used without a bath, and shown to be capable of producing Turing patterns after long simulation times. The significance of our work is that RMPC can serve as a feasible simulation tool for both short and long-term simulations, can handle spatial inhomogeneities, can model a fairly large range of diffusivities in a reaction-diffusion scenario, and is capable of producing Turing patterns. An advantage of this method includes more detailed system information in feasible simulation times.


    Host detection of viral pathogens by the innate immune response is an effective first line of defence against infection. In eukaryotic cells, the innate immune system relies heavily on protein recognition of evolutionarily conserved viral double-stranded RNA (dsRNA) to stimulate interferon production that prepares neighbouring cells for attack [1]. The interferon response then up-regulates the transcription of additional antiviral proteins to amplify the response. A key family of interferon-stimulated antiviral enzymes are the 2'-5'-oligoadenylate synthetases (OAS). OAS enzymes specifically bind viral dsRNA and through a complex cascade attenuate viral replication. An understanding of how these enzymes discriminate viral from host dsRNA remains a key unresolved question.

    OAS enzymes require viral dsRNA as a cofactor to polymerise ATP into chains of unusual 2'-5' linked oligoadenylates (2-5A) [2]. These chains bind to and activate RNase L, which degrades cellular and viral RNA thereby restricting viral propagation [3]. In humans, three genes (OAS1/2/3) each produce enzymes with one, two and three repeats of a 346 amino acid core OAS unit [4]. Each version has different dsRNA affinity, builds different 2-5A chain lengths, and responds to different viral infections. Interaction between multiple OAS enzymes (oligomerisation) appears to be central to these differences, presumably to respond to different viral dsRNA lengths. OAS1, comprising a single domain, forms monomer, dimers and tetramers, and will bind dsRNA of any length longer than 17 base pairs (bp) [5]. OAS2 can form a dimer and requires a minimum of 35 bp [6]. OAS3 is not known to oligomerize and requires at least 50 bp dsRNA for activation [7]. While high-resolution structures of OAS enzymes show that dsRNA is responsible for dramatic reorganisation of the enzyme active site, little is known regarding the impact of dsRNA binding on oligomerisation or the detailed catalytic mechanism.

    OAS enzymes are differentially expressed in specific cell types, and their presence in different sub-cellular fractions such as mitochondria, nuclear, and rough/smooth microsomal fractions suggest that these proteins may have distinct functional roles. OAS1, 2, and 3 differ significantly in their response to common members of the Flaviviridae family of human viral pathogens including Yellow fever, Dengue fever, Japanese encephalitis, Tick-borne encephalitis and West Nile virus. Preliminary data from [6] suggests that OAS1 and OAS2 may have distinct substrate specificity in terms of dsRNA length, which may explain their divergent roles in antiviral processes.

    The OAS enzymatic mechanism is complex as it involves dsRNA cofactors of variable length and structure, an ATP substrate that is polymerised into chains of variable length, and non-covalent OAS self-association. As a result of this complexity, the detailed biochemical process is not completely defined. In this work, we focus on OAS2 activation by dsRNA. Based on in vitro experimental data published in [6], the process appears to depend on concentrations and lengths of dsRNA; however, limitations of the experimental data that can only be accurately recorded at the endpoint of enzymatic activity prevents complete definition of the process. Hence, mathematical modelling is used to investigate possible OAS2 activation mechanisms and the effects of dsRNA lengths. Plausible biochemical scenarios are translated to mathematical models and their responses are compared to in vitro experimental data to test different hypotheses. In total, three kinetic models are derived. After calibration of models to data, a model selection method is used to determine the model that best represents the data. Furthermore, analysing parameter estimates allows us to propose a mechanism accounting for the role of dsRNA length in OAS2 activation.

    OAS2 enzymes require dsRNA as a cofactor for enzymatic activity. In the presence of substrate ATP, OAS2 binds to dsRNA and the enzyme polymerises ATP into chains of 2-5A leaving pyrophosphate (PPi) as a byproduct [2]. In vitro experimental data considered in this work are concentrations of PPi recorded at six time points (including the initial time, t0=0 min). The byproduct concentration is measured at t1=5 min, t2=10 min, t3=15 min, t4=20 min and t5=30 min for one concentration cO=0.3μg/mL of OAS2, and five concentrations of dsRNA are considered with seven different lengths [6]. The five dsRNA concentrations cj are c1=0.5μg/mL, c2=1μg/mL, c3=2μg/mL, c4=4μg/mL and c5=8μg/mL. The seven lengths Lk considered for dsRNA are L1=40bp, L2=50bp, L3=60bp, L4=70bp, L5=80bp, L6=90bp and L7=120bp. The ATP concentration used in the experiments is 1mg/mLcO. Experimental data used in this work have been published in [6]; for details on in vitro experiments yielding them refer to [6].

    In the in vitro experiments considered [6], the substrate, which is ATP, is in abundance. Hence, the contribution of substrate to enzyme kinetics is not of interest in this work; no equation describes ATP concentration in models. The OAS2 enzyme consists of two OAS domains [8]. As shown in Koul et al. [6], a minimal length of 35bp for dsRNA is required to activate OAS2, and both OAS2 domains must be bound simultaneously by dsRNA to lead to product formation. Furthermore, the binding and unbinding of enzymes and cofactors are reversible reactions; neither enzyme nor cofactor are consumed during these reactions nor the formation of product and byproduct. Finally, as previously mentioned, the experimental data only monitors the concentration over time of the byproduct (PPi) for a given enzyme concentration and different cofactor concentrations and lengths. In the experimental data, the formation of byproduct is sensitive to both cofactor concentrations and lengths.

    Based on these observations and available in vitro experimental data three mathematical models are developed to investigate the activation of OAS2 by dsRNA and the role of the dsRNA length on the activation (Figure 1).

    Figure 1.  Three models considered to investigate OAS2 activation by dsRNA.

    The first model considers an enzymatic kinetics describing only the kinetics of the enzyme E and cofactor R as the substrate is in abundance. The activated enzyme or productive stage D, which is formed by the OAS2 and the cofactor bound to the two domains of OAS2, is not consumed in the formation of byproduct P. Assumptions of Model-E are as follows (Figure 1)

    E+Rk1k1mDØkpDP

    and the model equations are

    dEdt=k1ER+k1mD,dRdt=k1ER+k1mD,dDdt=k1ERk1mD,dPdt=kpD, (3.1)

    where k1 is the rate of association of an OAS2 and dsRNA, it describes the rate of production of the compound D or activated enzyme. To reduce the model complexity, binding of the two domains of OAS2 to the same dsRNA is not explicitly described and is approximated as an one-step binding. The rate k1m is the dissociation rate of the enzyme and cofactor. Note that in Model-E once a dsRNA molecule is bound with an OAS2, it cannot accommodate another OAS2 enzyme. The parameter kp is the rate of formation of the byproduct by the activated enzyme. Parameters are listed in Table 1. Model-E is considered with non-negative initial conditions E(0)=cO, R(0)=cR, D(0)=0 and P(0)=0, where the positive constants cO and cR are the initial concentrations of OAS2 and dsRNA, respectively.

    Table 1.  Parameters used in models. (Top) Fixed parameters whose values are set from experimental conditions [6]. (Bottom) Free parameters whose values are estimated by calibrating models responses to experimental data.
    Parameter Description (units) Model-
    cO Initial concentration of OAS2 (0.3μg/mL) E, N, M
    cR Initial concentration of dsRNA (0.5, 1, 2, 4 and 8μg/mL) E, N, M
    L Length of dsRNA (40, 50, 60, 70, 80, 90 and 120 bp) M
    2RO Diameter of OAS2 (35 bp) M
    k1 Association rate of an OAS2 enzyme to an unoccupied dsRNA ((μ g/mL)1.min1) E, N, M
    k1m Dissociation rate of the OAS2 enzyme and dsRNA (min1) E, N, M
    ki Association rate of the ith enzyme to an dsRNA ((μ g/mL)1.min1) M with i={2,,n}, where n=L2RO
    kim Dissociation rate of ith enzyme from an dsRNA (min1) M with i={2,,n}, where n=L2RO
    k Association rate of a non-productive stage ((μ g/mL)1.min1) N
    km Dissociation rate of non-productive stage (min1) N
    kp Rate of byproduct formation (min1) E, N, M

     | Show Table
    DownLoad: CSV

    A way to represent the impact of the dsRNA length on OAS2 activation is to consider that some interactions are not complete [9] as the two domains of OAS2 have to bind to dsRNA to produce an active OAS2. Hence, the formation of compounds in which both domains are not properly attached to cofactor is hypothesised. These compounds are called non-productive compounds, their specific forms are not described. Based on this assumption, the non-productive binding stage N that cannot result in the formation of the product and byproduct is added to Model-E. After a non-productive binding, OAS2 can detach from dsRNA and attach to dsRNA again. The assumptions for Model-N are as follows (Figure 1)

    yielding the model equations

    dEdt=k1ER+k1mDkER+kmN,dRdt=k1ER+k1mDkER+kmN,dDdt=k1ERk1mD,dNdt=kERkmN,dPdt=kpD, (3.2)

    where the additional parameters k and km are the association and dissociation rates of the non-productive compound N, respectively (Table 1). Model-N is considered with the initial conditions E(0)=cO, R(0)=cR, D(0)=N(0)=0 and P(0)=0.

    Another possible effect of cofactor length on OAS2 activation is the occurrence of multi-binding of OAS2 on a same dsRNA. dsRNA is not consumed nor changed by the OAS2 activation; OAS2 interacts with dsRNA, facilitates product formation, and then detaches from dsRNA. In the in vitro experiments considered, the longer the dsRNA used, the more product formation observed at a given OAS2 concentration [6]. Hence, another reasonable scenario is that longer dsRNA able to accommodate more than one OAS2 (presumably > 70 bp) can be bound with more than one OAS2 simultaneously, leading to an increased product formation. Hence, the multi-binding assumption that allows up to n OAS2 bound to one dsRNA simultaneously, is summarised as follows (Figure 1):

    E+Rk1k1mD1+Ek2k2mD2+EDn1+EknknmDn+EkpD1PkpDiiPkpDnnP

    where n=L2RO, with L the length of dsRNA and RO the radius of an OAS2 enzyme. The compound Di represents i activated OAS2 attached to the same dsRNA with i{1,,n}. As shown in [6], the activity of OAS1, which is the monomeric equivalent of the catalytic domain of OAS2, does not depend on cofactor length. Therefore, the rate of formation of byproduct per activated enzyme is assumed to be constant and does not depend on the number of OAS2 bound to the same dsRNA. The compound Di produces i byproducts at rate kp.

    The kinetics of OAS2 in presence of dsRNA of length L is then described by the following system of n+3 equations:

    dEdt=k1ER+k1mD1+n1i=1[k(i+1)mDi+1k(i+1)DiE],dRdt=k1ER+k1mD1,dD1dt=k1ERk1mD1+k2mD2k2D1E,dDjdt=kjDj1EkjmDj+k(j+1)mDj+1k(j+1)DjE,j{2,3,,n1},dDndt=knDn1EknmDn,dPdt=kpni=1iDi, (3.3)

    where the additional parameters ki with i>1 are the rates of subsequent association of one OAS2 to a dsRNA already associated with i1 OAS2. The rates kim with i>1 are the dissociation rates of one OAS2 from a dsRNA bound with i OAS2. The initial conditions for Model-M are E(0)=cO, R(0)=cR, Di(0)=0 with i{1,,n} and P(0)=0. Note that when dsRNA used in experiments have a length L smaller or equal to 70bp, then Model-M is considered with n=1 which is equivalent to Model-E.

    Before calibrating models to data, their behaviour is characterised. Mathematical analysis of models is carried out and presented in [10]; all models are well-posed. Furthermore, reduced models obtained by decoupling the P equation from systems (3.1), (3.2) or (3.3) are shown to have an unique stable positive steady state (results not shown; see [10]).

    To calibrate models, numerical solutions of P in systems (3.1), (3.2) or (3.3) are compared to in vitro experimental data, published in [6], that track the concentration of PPi over time for different concentrations of dsRNA cj (=cR) of a given length Lk (=L). The experimental data are shown in Figures 2 or 5 as red errors bars.

    Figure 2.  Fitting results obtained with the first strategy (all concentrations for one length). (A) Model-E. (B) Model-N. (C) Model-M. In each panel, each subfigure shows the PPi concentrations (in μg/mL) over time (in minutes) for the five concentrations of dsRNA arranged in columns and the seven lengths arranged in rows. Red error bars and dots are the experimental data, the concentrations of PPi at the six time points. Grey curves represent all model solutions using parameter value estimates obtained with the global optimisation method. Blue curves are the best fits with errors RSSLk(ˆpBestˆpBest).
    Figure 3.  Parameter estimates with first strategy. Each dot represents one estimate ˆpˆp obtained by calibrating models with the heuristic optimisation method minimising (4.1). For Model-M and lengths from 80 to 120 bp, in the same graph, dots represent estimates in the k1m/k1kp plane and + in the k1m/k1k2m/k2 plane. Colours are defined as RSSLk(ˆpˆp)/max(RSSLk(ˆpˆp)). Blue codes for the smallest "relative error" (the most likely parameter values given the considered data) and red codes for largest values of relative errors. Units for kp are min1 and for dissociation constants km/k μg/mL.
    Figure 4.  Distributions of parameter values obtained with the first strategy of model calibration when each length is considered independently: (A) Model-E, (B) Model-N and (C) Model-M. Vertical lines in distributions indicate 95% confidence intervals for the given parameter. Units for kp are min1 and for dissociation constants km/k μg/mL.
    Figure 5.  Fitting results obtained with the second strategy (all concentrations for all length). (A) Model-E. (B) Model-N. (C) Model-M. In each panel, each subfigure shows the PPi concentrations (in μg/mL) over time (in minutes) for the five concentrations of dsRNA arranged in columns and the seven lengths arranged in rows. Red error bars and dots are the experimental data, the concentrations of PPi at the six time points. Grey curves represent all model solutions using parameter value estimates obtained with the global optimisation method. Blue curves are the best fits with errors RSS(ˆpBestˆpBest).

    Once each model is calibrated to experimental data, the best model amongst the three models considered is selected based on the data, goodness of the fit and model parsimony using the corrected Akaike Information Criterion (AICc) [11,12].

    Experimental data (in red in Figures 2 or 5) show that the formation of byproduct increases with the concentration and length of the cofactor. Although the aim of the work is to find a mechanism that can explain the formation of byproduct for all lengths and concentrations of dsRNA simultaneously, we wanted first to verify if models can accommodate changes in concentrations for a given length and if the effect of dsRNA length on the production of P can be identified through the variation of the parameter estimate values obtained when models are calibrated to each length independently. Thus, two strategies are used to estimate free parameter values listed in Table 1. First, model responses are fitted to experimental data obtained with the five concentrations of dsRNA for a given length at the five time points. The second strategy considers all five concentrations for seven lengths of dsRNA at the five time points.

    In the first strategy, for the length Lk with k{1,,7}, the cost function to be minimised is:

    RSSLk(pp)=mii=1(mjj=1(Pcj,Lk(pp,ti)Dti,cj,Lk)2), (4.1)

    where pp is the free parameter vector, Pcj,Lk(pp,ti) is the predicted concentration of PPi (solution for P of a given model) at time ti for the concentration of dsRNA cR=cj of length L=Lk. The experimental concentrations of PPi are Dti,cj,Lk at time ti for the concentration of dsRNA cR=cj of length L=Lk. The parameters mj(=5) and mi(=5) are the numbers of dsRNA concentrations and time points considered in the study, respectively. In this first calibration strategy, parameters are estimated on 25 observation points by minimising (4.1).

    In the second strategy, parameters are estimated over 175 observation points and the cost function is then:

    RSS(pp)=mkk=1(mii=1(mjj=1(Pcj,Lk(pp,ti)Dti,cj,Lk)2)), (4.2)

    where mk(=7) is the number of dsRNA length considered in the study.

    To find the estimates ˆpˆp for free parameter values for both strategies, the optimisation problem RSS(ˆpˆp)=minppRSS(pp) is numerically solved using the genetic algorithm developed in the R package GA [13,14]. To obtain confidence intervals for estimates ˆpˆp, the genetic algorithm is run several times and let converge to obtain several plausible points in the multidimensional parameter space. The likelihood of a particular estimate ˆpˆp is a measure of how well the data supports that particular value of parameters. Maximising the likelihood is equivalent to minimising the negative of the log-likelihood function. Assuming independent and normally distributed additive measurement errors, the negative log-likelihood of the estimates ˆpˆp for a given model can be roughly approximated by RSSLk(ˆpˆp) for the first or RSS(ˆpˆp) for the second calibration strategy [15]. Hence, the likelihood of different estimates obtained with the heuristic approach is assessed using RSS values.

    For each calibration strategy (one length or all lengths), the AICc of model i{ModelE,ModelN,ModelM} is computed as follows:

    AICci=Mln(RSSi(ˆpBestˆpBest)M)+2KiMMKi1, (4.3)

    where Ki is the number of estimated parameters (number of free parameters in mathematical model i + 1), M is the number of observations. For the one length strategy, M is equal to 25 whereas, for the all lengths strategy, M is 175. The constant RSSi(ˆpBestˆpBest) is the smallest value of the errors obtained in multiple runs of the genetic algorithm in calibrating model i for a given strategy. For each calibration strategy, the best model is the one with the smallest AICc.

    To help the interpretation of results, the Akaike weights are then computed [15,16]. The Akaike weight wi of model i is

    wi=exp(Δi2)Qq=1exp(Δq2), (4.4)

    where Δi=AICciminiAICci and Q(=3) is the number of models considered in the study. The Akaike weights can be interpreted as the probability of the models to be the most likely given the experimental data and set of models considered. A model can be considered as the single best model if its weight wi satisfies wi>0.95.

    Fitting results obtained with the first and second strategy are shown in Figures 2 and 5, respectively. Parameter estimates, which are obtained with the first strategy, are shown in multidimensional parameter spaces with information related to their likelihood in Figure 3 and distributions and confidence intervals (95% of values obtained by calibrating models) of parameter values are shown in Figure 4. Details on parameter estimates obtained with the second calibration strategy are given in Figure 6. Results of model selection for both strategies are given in Table 2.

    Figure 6.  Results of model calibration to all concentrations and lengths data: (first row) Model-E, (second row) Model-N and (third row) Model-M. (First and second columns) Parameter estimates in multidimensional parameter spaces: each dot represents one of estimates ˆpˆp obtained by minimising (4.2). Colors, which are defined as RSS(ˆpˆp)/max(RSS(ˆpˆp)), are used to assess estimates likelihood. Blue codes for the smallest "relative error" (the most likely parameter values given the considered data) and red codes for largest values of relative errors. (Last column) Distributions of parameter values with 95% confidence intervals delimited by vertical lines. Units for kp are min1 and for dissociation constants km/k μg/mL.
    Table 2.  Model selection results. RSS(ˆpBestˆpBest) are the smallest errors obtained among the several runs of the genetic algorithm to minimise (4.1) for the one length calibration strategy or (4.2) for the all lengths calibration strategy. K is the number of estimated parameters used in the computation of AICc defined in (4.3) (number of parameter in the mathematical model + 1). AICc weights for each model and calibration strategy are computed using (4.4). Bold indicates the best models.
    dsRNA length 40bp 50bp 60bp 70bp 80bp 90bp 120bp All
    RSS(ˆpBest) Model-E 6.133 55.40 86.56 198.3 330.5 372.0 276.8 30753.3
    Model-N 6.113 55.27 86.56 200.6 352.1 375.4 341.0 30779.6
    Model-M 6.131 55.40 86.56 199.9 247.6 341.6 234.7 9685.4
    K Model-E 4 4 4 4 4 4 4 4
    Model-N 6 6 6 6 6 6 6 6
    Model-M 4 4 4 4 6 6 8 8
    AICc Model-E 25.13 29.89 41.05 61.77 74.54 77.50 70.11 912.80
    Model-N 18.54 36.50 47.71 68.73 82.79 84.39 81.99 917.22
    Model-M 25.14 29.89 41.05 61.98 73.99 82.03 80.98 719.25
    AICc weight Model-E 0.49 0.49 0.49 0.52 0.43 0.88 0.993 9×1043
    Model-N 0.02 0.02 0.02 0.01 0.01 0.03 0.003 1×1043
    Model-M 0.49 0.49 0.49 0.47 0.56 0.09 0.004 1

     | Show Table
    DownLoad: CSV

    When each length of dsRNA is considered independently, all models well represent the formation of PPi over time for all concentrations of dsRNA of a given length (Figure 2). Models accommodate changes in the cofactor concentration.

    For the three models, values of the rate of production of PPi, kp, are getting larger as the dsRNA length increases (Figures 3 and 4). When only one OAS2 is allowed on a single dsRNA molecule, there exists a linear relationship between the dissociation constant of enzyme and dsRNA k1m/k1 and the rate kp (Figure 3). Interestingly, for the shortest length 40bp, values for k1m/k1 are larger than kp values whereas this trend is reversed for lengths longer than 60bp. For short dsRNA, the binding of enzyme and cofactor is the limiting mechanism whereas for the longer dsRNA the limiting factor is the formation of product.

    Model-E appears to be the most identifiable. Intervals of parameter estimates obtained for Model-E are quite narrow and fitting errors do not vary much (Figures 3 and 4A). The genetic algorithm always converges to similar parameter values with similar errors when calibrating Model-E for individual length data.

    On the contrary, intervals for parameter estimates obtained for Model-N are more disperse and errors span over wider ranges (Figures 3 and Figure 4B). For Model-N, similar parameter value distributions (Figure 4B) with more clustered estimates in parameter spaces (Figure 3) are obtained as the length gets larger than 70bp. Interestingly, for long lengths (more than 80bp), the formation of non-productive compounds is mostly found less favourable than the formation of productive compounds (km/k>k1m/k1) (Figure 3). For long lengths (more than 80bp), unproductive states are short-lived in comparison to productive states.

    For Model-M, when dsRNA are long enough to allow the binding of multiple OAS2 to the same dsRNA (from 80 to 120bp), we can observe that the dissociation constants kim/ki for the subsequent bindings of OAS2 are smaller than those of prior ones, k1m/k1>k2m/k2>k3m/k3 (Figures 3 and 4C). This result indicates that subsequent bindings are more favourable than the prior ones, which suggests a cooperative dynamics for OAS2 activation.

    Keeping in mind that for lengths smaller or equal to 70bp Model-M and Model-E are equivalent, Model-E is found to be the best model using AICc when dsRNA are shorter than or equal to 70bp (Table 2). For 80 and 90bp long dsRNA, the confidence set of models is Model-E and Model-M. Finally, for the 120bp length, the simplest model Model-E is selected as the unique best model. However, note that if only the goodness of fit (RSS(ˆpBestˆpBest)=RSSLk(ˆpBestˆpBest) in Table 2) is used as a criterion the most complex model Model-M is the best model for dsRNA longer than 70bp, whereas the three models are indistinguishable for dsRNA shorter or equal to 70bp.

    When all dsRNA lengths are considered simultaneously, Model-E and Model-N well represent the experimentally observed formation of PPi over time for all concentrations of dsRNA of intermediate lengths 60-80bp, whereas they fail to accurately represent byproduct formation for shorter and longer lengths (Figure 5A-B). As with the first calibration strategy, estimates for shared dissociation constants and kp are found to be of the same order of magnitude for Model-E and Model-N (Figure 6). Despite the additional experimental data considered with the second calibration strategy, Model-N still appears to be highly non-identifiable.

    On the other hand, Model-M gives a good representation of all dsRNA concentrations regardless of dsRNA length (Figure 5C). Despite requiring the largest number of parameters, Model-M is also selected as the unique best model when compared to all data (Table 2). Furthermore, note that the use of additional data in the second strategy improves the identifiability of Model-M: when considering all experimental data, narrower confidence intervals for estimates are obtained and estimate values are more clustered in multidimensional parameter spaces (Figure 6).

    Finally, for the best model, Model-M, all estimates found for k1m/k1 and kp are linked by the same linear relationship with mostly kp<k1m/k1. Furthermore, there exists the same linear relationship between all estimates found for k1m/k1 and k2m/k2 ensuring k1m/k1>k2m/k2 (Figure 6). Estimates for dissociation constants of the activation by the same dsRNA of a second (k2m/k2) and a third (k3m/k3) OAS2 have similar values. The dissociation constant of the first enzyme (k1m/k1) is more than ten times larger than the second (k2m/k2) and third (k3m/k3). The second and third enzymes have higher affinity (which is the inverse of the dissociation constant) than the first enzyme. This suggests that subsequent bindings of OAS2 are facilitated by the first binding and that there is cooperative binding of OAS2 to a single dsRNA molecule.

    In this work, using mathematical modelling we investigate the effect of dsRNA length on the activation of OAS2; dsRNA serve as cofactors to enable the enzymatic activity of OAS2. Based on experimental data, our working hypothesis is that the activity of OAS2 increases as its cofactor length increases. Three models are designed: a simple enzymatic activity (Model-E), enzymatic activity with reversible inhibition (Model-N) and activation of multiple OAS2 by a single dsRNA (Model-M). In Model-E, the cofactor length is not explicitly accounted for, whereas Model-N considers the possibility of the formation of a reversible unproductive RNA-protein interaction that results from the binding of only one domain of OAS2 to dsRNA due to its length. Model-M more explicitly incorporates the cofactor length by allowing a single dsRNA of sufficient length (more than 70bp) to bind to and activate multiple OAS2 enzymes. Moreover, two strategies are adopted to calibrate models. First, dsRNA lengths data are considered independently that yields an implicit length-dependency of the model parameters and gives us another tool to characterise the effect of the cofactor length on the different parts of the kinetics. In the second calibration strategy, all dsRNA concentrations and lengths data are considered simultaneously as the aim of our study is to propose a scenario that accommodates changes in both cofactor concentrations and lengths.

    The power/strength of our study lies on the consideration of multiple models tested against experimental data with various strategies and then compared via model selection; see also [17,18]. When considering kinetics for a single dsRNA length, the simplest enzymatic model Model-E is selected as the best model; however, when considering the kinetics for all dsRNA lengths simultaneously, the multiple binding model Model-M is the best model. Interestingly, a closer look at the estimates of parameters found for Model-M suggests that the multiple binding of OAS2 to the same dsRNA follows a cooperative regime (affinity of the activation of subsequent OAS2 enzymes is much higher than this of the first activation). Model-N, with transient inhibition, is never found to best explain the data even though it gives a good representation when lengths are considered independently. Note that a more complex version of Model-N with transitions from the non-productive compound N to productive compound D was also considered. As it provided similar results than Model-N but was even more penalized by the extra parameters, we did not present it in the work. However, the consideration of the formation of unproductive stages and its poor performance allow us to conclude that this scenario is not a mechanism that drives OAS2 activation by dsRNA.

    Hence, based on our analysis, we can propose a mechanism for OAS2 activation that accounts for the impact of the length of its cofactor: multiple OAS2 can bind to the same dsRNA, when its length permits, and subsequent bindings are facilitated which underlies an enzyme with a cofactor cooperativity. This hypothesis now has to be experimentally verified, but this seems to be in preliminary agreement with the recent observation emphasising stable self-association of OAS2 as central to its mechanism [8].

    Here, only three models are presented but more models have been developed in the base for this paper [10]. To the best of our knowledge, these mathematical models are the first models developed for the activation of OAS2 by dsRNA. Finally, other modelling approaches can be considered to account for the effect of cofactor length on the enzyme activation such as the use of length-dependent binding rate for enzymes to cofactors [19,20].

    The work was originally part of DL's MSc thesis at The University of Manitoba. SMK and SP are supported in part by NSERC Discovery Grants and a Faculty of Science Collaborative Grant (University of Manitoba).

    All authors declare no conflicts of interest in this paper.



    [1] D. T. Gillespie, Exact stochastic simulation of coupled chemical reactions, J. Phys. Chem., 81 (1977), 2340-2361. doi: 10.1021/j100540a008
    [2] D. T. Gillespie, Approximate accelerated stochastic simulation of reacting systems, J. Phys. Chem., 115 (2001), 1716-1733. doi: 10.1063/1.1378322
    [3] A. Stundzia, C. Lumsden, Stochastic simulation of coupled reaction-diffusion processes, J. Comput. Phys., 127 (1996), 196-207. doi: 10.1006/jcph.1996.0168
    [4] S. Isaacson, C. Peskin, Incorporating diffusion in complex geometries into stochastic chemical kinetics simulations, SIAM J. Sci. Comput., 28 (2006), 47-74. doi: 10.1137/040605060
    [5] D. Fange, J. Elf, Noise-induced min phenotypes in E. coli, PLoS Comput. Biol., 2 (2006), e80.
    [6] L. Ferm, A. Hellander, P. Lotstedt, An adaptive algorithm for simulation of stochastic reaction-diffusion processes, J. Comput. Phys., 229 (2010), 343-360. doi: 10.1016/j.jcp.2009.09.030
    [7] J. M. A. Padgett, S. Ilie, An adaptive tau-leaping method for stochastic simulations of reaction-diffusion systems, AIP Adv., 6 (2016), 035217.
    [8] R. Strehl, S. Ilie, An adaptive tau-leaping method for stochastic systems with slow and fast dynamics, J. Chem. Phys., 143 (2015), 234108.
    [9] C. A. Smith, C. A. Yates, Spatially extended hybrid methods: A review, J. R. Soc. Interface, 15 (2018), 20170931.
    [10] M. B. Flegg, S. J. Chapman, R. Erban, The two regime method for optimizing stochastic reaction-diffusion simulations, J. R. Soc. Interface, 9 (2012), 859-868. doi: 10.1098/rsif.2011.0574
    [11] A. Hellander, S. Hellander, P. Lotstedt, Coupled mesoscopic and microscopic simulation of stochastic reaction-diffusion processes in mixed dimensions, Multiscale Model. Sim., 10 (2012), 585-611. doi: 10.1137/110832148
    [12] M. B. Flegg, S. Hellander, R. Erban, Convergence of methods for coupling of microscopic and mesoscopic reaction-diffusion simulations, J. Chem. Phys., 289 (2015), 1-17.
    [13] J. Hattne, D. Fange, J. Elf, Stochastic reaction-diffusion simulation with MesoRD, Bioinformatics, 21 (2005), 2923-2924. doi: 10.1093/bioinformatics/bti431
    [14] S. S. Andrews, D. Bray, Stochastic simulation of chemical reactions with spatial resolution and single molecule detail, Phys. Biol., 1 (2004), 137-151. doi: 10.1088/1478-3967/1/3/001
    [15] J. S. van Zon, P. R. ten Wolde, Greens-function reaction dynamics: A particle-based approach for simulating biochemical networks in time and space, J. Chem. Phys., 123 (2005), 234910.
    [16] J. S. van Zon, P. R. ten Wolde, Simulating biochemical networks at the particle level and in time and space: Green's function reaction dynamics, Phys. Rev. Lett., 94 (2005), 128103.
    [17] S. J. Chapman, R. Erban, S. A. Isaacson, Reactive boundary conditions as limits of interaction potentials for Brownian and Langevin dynamics, SIAM J. Appl. Math., 76 (2016), 368-390. doi: 10.1137/15M1030662
    [18] H. G. Othmer, S. R. Dunbar, W. Alt, Models of dispersal in biological systems, J. Math. Biol., 26 (1988), 263-298. doi: 10.1007/BF00277392
    [19] Y. Cao, R. Erban, Stochastic Turing patterns: Analysis of compartment-based approaches, B. Math. Biol., 76 (2014), 3051-3069. doi: 10.1007/s11538-014-0044-6
    [20] M. B. Flegg, Smoluchowski reaction kinetics for reactions of any order, SIAM J. Appl. Math., 76 (2016), 1403-1432. doi: 10.1137/15M1030509
    [21] A. Malevanets, R. Kapral, Mesoscopic model for solvent dynamics, J. Chem. Phys., 110 (1999), 8605-8613. doi: 10.1063/1.478857
    [22] A. Malevanets, R. Kapral, Solute molecular dynamics in a mesoscale solvent, J. Chem. Phys., 112 (2000), 7260-7269. doi: 10.1063/1.481289
    [23] K. Tucci, R. Kapral, Mesoscopic model for diffusion influenced reaction dynamics, J. Chem. Phys., 120 (2004), 8262-8270. doi: 10.1063/1.1690244
    [24] K. Tucci, R. Kapral, Mesoscopic multiparticle collision dynamics of reaction-diffusion fronts, J. Chem. Phys. B, 109 (2005), 21300-21304. doi: 10.1021/jp052701u
    [25] K. Rohlf, S. Fraser, R. Kapral, Reactive multiparticle collision dynamics, Comput. Phys. Commun., 179 (2008), 132-139. doi: 10.1016/j.cpc.2008.01.027
    [26] K. Rohlf, Stochastic phase-space description for reactions that change particle numbers, J. Math. Chem., 45 (2009), 141-160. doi: 10.1007/s10910-008-9373-8
    [27] J. M. Yeomans, Mesoscale simulations: Lattice Boltzmann and particle algorithms, Physica A, 369 (2006), 159.
    [28] P. J. Hoogerbrugge, J. M. V. A. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, EPL, 19 (1992), 155.
    [29] S. Bedkihal, J. C. Kumaradas, K. Rohlf, Steady flow through a constricted cylinder by multiparticle collision dynamics, Biomech. Model. Mechan., 12 (2013), 929-939. doi: 10.1007/s10237-012-0454-z
    [30] T. Akhter, K. Rohlf, Quantifying compressibility and slip in multiparticle collision (MPC) flow through a local constriction, Entropy, 16 (2014), 418-442. doi: 10.3390/e16010418
    [31] K. Rohlf, Compressible slip flow through constricted cylinders with density-dependent viscosity, Dynam. Cont. Dis. Ser. B, 25 (2018), 233-257.
    [32] T. Ihle, D. M. Kroll, Stochastic rotation dynamics: A Galilean-invariant mesoscopic model for fluid flow, Phys. Rev. E, 63 (2001), 020201.
    [33] T. Ihle, D. M. Kroll, Stochastic rotation dynamics: I. Formalism, Galilean invariance, and Green-Kubo relations, Phys. Rev. E, 67 (2003), 066705.
    [34] T. Ihle, D. M. Kroll, Stochastic rotation dynamics: II. Transport coefficients, numerics, and long-time tails, Phys. Rev. E, 67 (2003), 066706.
    [35] H. Noguchi, G. Gompper, Transport coefficients of off-lattice mesoscale-hydrodynamics simulation techniques, Phys. Rev. E, 78 (2008), 016706.
    [36] E. Tüzel, T. Ihle, D. M. Kroll, Dynamic correlations in stochastic rotation dynamics, Phys. Rev. E, 74 (2006), 056702.
    [37] R. Strehl, K. Rohlf, Multiparticle collision dynamics for diffusion-influenced signaling pathways, Phys. Biol., 13 (2016), 046004.
    [38] A. Sayyidmousavi, K. Rohlf, Reactive multi-particle collision dynamics with reactive boundary conditions, Phys. Biol., 15 (2018), 046007.
    [39] V. Mortazavi, M. Nosonovsky, Friction-induced pattern formation and Turing systems, Langmuir, 27 (2011), 4772-4779. doi: 10.1021/la200272x
    [40] D. A. Garzon-Alvarado, C. H. Galeano, J. M. Mantilla, Computational examples of reaction-convection-diffusion equations solution under the influence of fluid flow: A first example, Appl. Math. Model., 36 (2012), 5029-5045. doi: 10.1016/j.apm.2011.12.041
    [41] J. Wei, M. Winter, Flow-distributed spikes for Schnakenberg kinetics, J. Math. Biol., 64 (2012), 211-254. doi: 10.1007/s00285-011-0412-x
  • This article has been cited by:

    1. Deokro Lee, Amit Koul, Nikhat Lubna, Sean A. McKenna, Stéphanie Portet, Mathematical modelling of OAS2 activation by dsRNA and effects of dsRNA lengths, 2021, 6, 2473-6988, 5924, 10.3934/math.2021351
    2. Victor Ogesa Juma, Leif Dehmelt, Stéphanie Portet, Anotida Madzvamuse, A mathematical analysis of an activator-inhibitor Rho GTPase model, 2022, 9, 2158-2491, 133, 10.3934/jcd.2021024
    3. Marc R. Roussel, Moisés Santillán, Biochemical Problems, Mathematical Solutions, 2022, 7, 2473-6988, 5662, 10.3934/math.2022313
  • Reader Comments
  • © 2019 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(4923) PDF downloads(351) Cited by(1)

Figures and Tables

Figures(15)  /  Tables(2)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog