BACKGROUND
computer models contribute to the integrative understanding of complex molecular processes in the cell.
the most commonly used approach is deterministic modeling based on ordinary differential equations .
when the studied system comprises species with a low molecular count, stochastic formalisms, e.g.
gillespie’s algorithm that simulates trajectories and uses discrete molecule numbers are more appropriate  <cit> .
in recent years, the use of stochastic models has substantially increased  <cit> .
considering stochasticity in biological systems has changed the quantitative and qualitative understanding obtained by previous deterministic models  <cit> .
examples of biological phenomena discovered by stochastic modeling include gene expression in burst-like patterns  <cit> , productive or latent cell decision after hiv-infection  <cit> , and the presence of oscillatory behavior induced by noise  <cit> .
however, in contrast to the plethora of methods used to fit and analyze deterministic models , there are only very limited sets of methods available to do the same with stochastic models.
thus, only recently, methods for parameter fitting of stochastic models have been developed and are still under development.
parameter estimation methods for stochastic models have been so far designed for single time course data .
for experiments leading to data distributions, they include the moment closure , and the comparison of experimental and simulation distributions .
despite the current efforts, there are still major problems that limit the full applicability of those methods.
the main drawback is the high computational execution time.
to solve this problem novel efficient strategies to fit stochastic models are needed.
in the following we introduce a new approach that is a variant of existing methods based on the comparison of distributions.
differently to the existing methods we use the experimental data to define a condition that the model must fulfill in the deterministic regime.
this condition tests, if the deterministic steady states are in the close proximity to the modes of the experimental distribution.
only if the model fulfills this condition, the parameter set is evaluated in the stochastic regime.
in this way, the algorithm directs the evaluation of parameter sets towards regions in parameter space with high potential to reproduce the experimental data.
additionally, this method applies to non-linear models and complex experimental data distributions.
to test the potential of this new method, we selected an open scientific question and real experimental data.
we investigated the molecular mechanisms responsible for the experimentally observed bimodal dynamics of irf <dig> expression after interferon  stimulation.
it is well documented that in isogenic cell populations not all cells produce an antiviral response when infected by a virus  <cit> .
the mechanisms behind this heterogeneity have been related to stochastic events in the signaling pathways responsible to elicit the antiviral response .
ifn is a cytokine that activates the jak-stat signaling pathway in virus-infected cells.
the jak-stat signaling pathway induces the translocation of isgf <dig>  into the nucleus to directly activate the transcription of a set of several hundred ifn-stimulated genes   <cit> .
irf <dig> is an isg with a central role in the immune response  <cit> .
recently, it was observed that after ifn stimulation murine fibroblasts show a switch-like pattern of irf <dig> expression, which is reflected at the population level as bimodality  <cit> .
we developed a model to describe the observed bimodal dynamics of irf <dig> expression after ifn stimulation.
we hypothesize that the binding of irf <dig> and isgf <dig> to the irf <dig> promoter is the mechanism responsible for this bimodal behavior.
our simulation results reproduced irf <dig> switch-like dynamics at single cell level, and bimodality was achieved at the population level.
hereby, we used the newly developed method to fit the model and correctly reproduce the experimental data.
our results allow us to conclude that the irf <dig> promoter activation by the combination of irf <dig> dimers and isgf <dig> is sufficient to quantitatively explain the observed irf <dig> bimodal dynamics.
methods
experimental data
published experimental data were produced by rand et al.
<cit> .
the data described the expression of irf <dig> after ifn- β stimulation in a population of murine nih3t <dig> fibroblasts.
rand’s experiments were done in the following way: first, cells were transfected with a bac  containing irf <dig> and reporter mcherry genes fused, subsequently cultures were treated with different concentrations of murine ifn- β. for illustrative purposes we selected the treatment with 150u of ifn- β, a concentration where bimodality is prevalent.
then, fluorescence in the cultures was monitored using flow cytometry at different time points during  <dig> hours.
numerical methods
in the deterministic regime the model was simulated using symbolic methods in matlab  <cit>  and/or using the lsoda algorithm in copasi  <dig>   <cit> .
stability was calculated by making the right-hand side of the differential equations equal to zero and subsequently by finding the roots of the system using function solve in matlab.
eigenvalues were calculated using function eig in matlab.
for solving the model under stochastic dynamics we used the gibson and bruck algorithm  <cit>  coded in copasi.
the random search and genetic algorithm were coded in matlab.
raw experimental data was analyzed using the function fcs data reader coded in matlab  <cit> .
the modes in the distributions were calculated using the function peakfinder coded in matlab  <cit> .
the source code of the project can be accessed via: https://sourceforge.net/projects/irf7-bimodaldynamics/.
RESULTS
new algorithm to fit stochastic biological models
comparing experimental and simulation distributions
the measurements of fluorescence were made comparable with the corresponding observable chemical species s
o in the model by a function h that maps state s
o to observation s
† as follows:
<dig> s†=h),
at each measurement time point t
i, i= <dig> …,n.
subsequently, a process to compare distributions was developed based on  <cit> .
first, considering a specific set of parameter values θ={θ
<dig> …,θ
d}, we performed ns repetitions of the stochastic simulations s={s
<dig> …,s
ns}.
the total of those stochastic simulations were used to build histograms with a fixed number of bins l for each t
i. subsequently, the  probability density function  of the simulations p
s,θ,b
l) were computed using the center of each bin b
l, for l= <dig> …,l, and the normalized form of the histograms.
on the other hand, having nm repetitions of single cell experimental data m={m
<dig> …,m
nm} histograms and pdfs p
e,b
l) were built.
even though, complex formulas can be use to calculate the number of bins in an histogram, here l was kept constant for p
s and p
e and was calculated as l=nm, this is a simple approach used by default by most programming languages .
finally, an approach using the difference of squares was used as a metric to calculate the distance between p
s,θ,b
l) and p
e,b
l), as follows:
<dig> f=∑i=1n∑l=1l,bl)−ps,θ,bl)) <dig>
algorithm implementation
the objective function f can be used to calculate a parameter estimate
<dig> θ^=argminθf.
as optimization usually requires plenty of function evaluations and each evaluation of f requires ns stochastic simulations, a direct optimization strategy is computationally unfeasible even for simple models.
therefore, we introduce a new strategy that selects good candidate parameters and only performs the stochastic simulations for these parameters.
in large systems where fluctuations can be discarded, the stochastic system can be reduced to the deterministic one  <cit> .
for this reason, in most cases the deterministic dynamics can be associated with a measure of central tendency in the pdfs obtained after the stochastic simulations.
using this reasoning we will introduce the strategy to efficiently estimate parameters for stochastic models making use of deterministic dynamics as an initial indicator.
assuming that the experimental data is in equilibrium at the last measurement point t
n, we used the modes from p
e,b
l) as a central measure of tendency, denoting the modes by α=,…,α
q), with q being the total number of modes.
then, using the ode version of the model ) we determine its stability.
here we tested whether the system has stable steady states and denote them with x∗=x1∗,…,xss∗, being s
s the number of stable steady states.
if the system has no stable steady state, it holds s
s= <dig> and the vector x
∗ is empty.
the calculation of steady states can either be performed by solving y˙= <dig> where y˙ is the right hand side of the ode system with x replaced by y. solving this equation can be impossible for large systems.
alternatively, one can numerically solve the ode systems until the flux is zero for all components.
varying the initial conditions leads to the different steady states.
the stability of the system was calculated after determining the sign of the eigenvalues , for stable steady states  <cit> .
having x
∗ and α we introduce the deterministic precondition:
<dig> βklowxk∗≤αk≤βkupxk∗,for1≤k≤q,and0≤βklow≤1and1≤βkup.
equation  means that the number of modes equals the number of stable steady states: s
s=q, and that each of the stable steady states is close to its corresponding mode.
only, if this is fulfilled, we calculate stochastic simulations and evaluate the objective function f. if the precondition is not fulfilled, we do not carry out stochastic simulations.
in this case, we assign a high  objective function value to direct the parameter search towards parameters that pass the deterministic precondition.
a graphical representation of the deterministic precondition is given in fig.
<dig>
fig.
<dig> the deterministic precondition.
left side the deterministic dynamics of the system can be seen, here the observable variable evolves to a single stable steady state x1∗.
this steady state can be calculated by finding the point in the system where the net flux is equal to zero x˙= <dig>  at the right hand side of the plot, the pdf of the experimental data is shown.
this pdf is rotated and the x-axis shows the count and the y-axis shows the bin values with units of molecules/cell.
from both graphs it can be seen that the deterministic precondition is valid when x1∗ is inside the range defined by β1lowx1∗ and β1upx1∗
at this point, it is important to notice that the deterministic precondition is not stating that the stable steady states form the ode system are equal to the mean obtained by solving the stochastic system.
rather, the deterministic precondition tests if the deterministic steady state lies around a certain range.
<dig> fcond=f,βklowxk∗≤αk≤βkupxk∗−−′′condition passed′′∞,else –′′condition not passed′′
our extended objective function with the deterministic check reads as:
and the estimate is defined as:
<dig> θ^cond=argminθfcond.
as previously discussed, our method is based on the central assumption that the modes in the experimental data are related to the stable steady states in the deterministic mathematical model.
if this assumption should fail, either of the following could happen: a) the method accepts a parameter that passed the deterministic precondition but does not show a good agreement of the distributions.
if it does not show good agreement, the objective function f will show a high value which means that this parameter does not lead to a good fit.
therefore, this case leads to a loss of computational time, but not to wrong results and it overall is still less expensive than evaluating every parameter set.
b) we reject a parameter that fails the deterministic precondition but would have lead to a good fit.
in this case, we lose indeed a good parameter.
this shows that our method is obviously a heuristic strategy that does not guarantee to find the exact global minimum, but this is true for almost any optimization strategy.
the algorithmic steps are depicted in fig.
<dig>
fig.
<dig> algorithm to fit stochastic models to experimental data.
the algorithm solves the model in the deterministic and stochastic regimes.
a condition observed in the experimental data is defined.
a set of parameter values is evaluated in the deterministic regime to test if the model reproduces this condition.
if the condition is met the stochastic simulation is performed.
otherwise, the parameter values are rejected.
the pdf obtained from the experimental data is compared with the pdf obtained after running the stochastic simulations.
this comparison is made using a difference of squares as an objective function.
this process is repeated until evaluating a total number of parameter sets or after a termination criterion is met.
the parameter set that best reproduces the experimental data is given by the minimum value obtained after the iterative evaluation of the objective function
parameter estimation methods
the deterministic precondition was first implemented using a random search algorithm.
random search is a global optimization strategy that tests random combinations of parameter values.
the successful output in this method is dependent on the total number of evaluated parameter sets  <cit> .
the more evaluated parameter values, the higher the probability to find the global minimum.
the pseudo-code for the random search is given in algorithm  <dig>
the second implementation of the deterministic precondition was using the more directed genetic algorithm.
genetic algorithm mimics evolution and is based on reproduction and selection.
this algorithm is made of a population of individuals , and each contains a genome that is defined by the number of parameters to optimize.
the individuals are ranked after solving the objective function, and a population of parental individuals is selected according to an elitism rate .
new individuals  are generated by pairing and recombining the parental genomes .
variability is introduced in the population by adding mutations in the new individuals according to a given mutation rate .
by the continuous process of selecting the best parameters after each generation, the algorithm evolves towards the regions in the parameter space that reduce the values in the objective function.
the pseudo-code for the genetic algorithm is given in additional file  <dig>
identifying parameters for a constitutive gene expression circuit with in-silico data
to prove the functionality and benefits of the new proposed algorithm, we first applied it to a simple example where the parameters are known a priori.
the model describes the stochastic dynamics of two variables, protein and mrna of a gene with constitutive expression  <cit> , a graphical representation of the constitutive gene expression circuit is given in fig.
3-a.
the model is described by the following reactions:
fig.
<dig> parameter estimation using a constitutive gene expression circuit and in-silico data.
a diagram of a gene with constitutive expression.
the mrna is produced at constant rate θ
<dig> and is degraded at constant rate θ
<dig>  the proteins are produced at constant rate θ
<dig> and are degraded with θ
<dig>  b comparison between in-silico data and model dynamics using the true parameters θ
o= and the estimated parameters θ^=.
a distribution with the in-silico data is given in grey.
in red is given the model dynamics.
c comparison between a priori and a posteriori parameter distributions.
the true parameters are represented by the intersection of the red lines.
in the priori distribution it can be observed that  <dig> parameters are randomly distributed in the parameter space.
the posteriori distributions were calculated by running  <dig> independent random searches with  <dig> parameters each, and by plotting the best parameter value selected by the algorithm.
in the posteriori distribution it is plotted only the parameters that passed the deterministic precondition, and it can be observed that those parameter estimates are in a region close to the true parameter values
<dig> ∅→θ1mrna,
<dig> mrna→θ2∅,
<dig> mrna→θ3protein,
<dig> protein→θ4∅.
to estimate the parameter values for this model we carried out the following procedure: first, we generated the in-silico data by selecting the following true parameter values θ
o= and running  <dig> independent stochastic simulations from which pdfs were built.
four observable time points were defined at  <dig> min,  <dig> min,  <dig> min and  <dig> min .
from the pdfs a mode was obtained at α
<dig> =  <cit>  a.u., where t
n= <dig> min.
then, we defined the deterministic precondition using eq.
, and assigning minimum and maximum acceptance ranges by setting β1low =  <dig>  and β1up =  <dig> , respectively.
a deeper analysis of the stability of the constitutive gene expression circuit is given in additional file  <dig>  to build simulated pdfs we used  <dig> repetitions of the stochastic model, this number was empirically calculated according to additional file  <dig>  for illustrative purposes, we assumed that the values for the parameters responsible for the mrna transcription and protein translation  were unknown, and the new algorithm was used to estimate those parameters.
the deterministic precondition was applied using the rs strategy obtaining that the algorithm only evaluates stochastic dynamics in  <dig> % of the tested parameter values, reducing in this way the total simulation time.
additionally, the complete algorithm was repeated  <dig> times and histograms of the estimated parameter were computed to determine whether they are close to θ
.
as can be observed in fig.
3-c the deterministic precondition reduces the evaluation of different parameters under stochastic dynamics by selecting only those parameters that are in a well-defined region in the proximity of θ
.
for this model, the main benefits of using the deterministic precondition was the reduction of the number of parameters evaluated under stochastic dynamics.
this rejection of parameters was made in an area outside the true parameter values, and hence no difference in accuracy is expected in comparison with a method without using the deterministic precondition.
a complete description of the analysis of the performance, accuracy and error for this example is given in additional file  <dig>  the model for the constitutive gene expression circuit can be obtained from biomodels database under reference model <dig>
mathematical model for irf <dig> expression dynamics
subsequently we applied our algorithm to a real problem with flow cytometry data.
here we studied the dynamics of murine irf <dig> gene expression upon ifn stimulation.
for this reason, we developed a model that comprises known key components and feedback mechanisms.
the overall system describes the active irf <dig> promoter  by the binding of irf <dig> dimer and isgf <dig> to the dna binding sites isre and irfe, respectively  <cit> , the transcription and translation of irf <dig>  and its subsequent phosphorylation and dimerization.
irf <dig> protein binding to the irfe binding site in the promoter results in the production of more irf <dig> protein, constituting a positive feedback loop  <cit> .
a graphical representation of the irf <dig> gene expression dynamics is given in fig.
<dig>
fig.
<dig> irf <dig> gene expression circuit.
a in viral infected cells irf <dig> gene expression is induced after ifn stimulation by the jak-stat signaling pathway.
the irf <dig> promoter activation is governed by the binding of irf <dig> dimer and isgf <dig> to the promoter dna binding sites .
irf <dig> promoter activation leads to the transcription of irf <dig> mrna and subsequent its translation to produce the irf <dig> protein.
notice that the irf <dig> production and subsequently phosphorylation and dimerization leads to the binding to its own dna binding site, resulting in the production of more irf <dig>  making in this way a positive feedback loop.
b a simplified system representing the irf <dig> gene expression dynamics is given.
here, the promoter transitions between active/basal states , the gene transcription and translation processes are represented as solid black lines as well as the feedback loop in the system
chemical reactions
the developed model consists of  <dig> species and  <dig> reactions  to ).
in the model, reaction  describes the ifn activation of the jak-stat signaling pathway.
for the description of this reaction and according to  <cit>  a saturable function was used.
downstream reactions of the pathway were lumped in reaction , so that the same variable was used to describe the output of the jak-stat signaling pathway, namely isgf <dig>  in this highly simplified model.
reaction  describes the irf <dig> promoter activation by the binding of irf <dig> dimer and isgf <dig> to isre and irfe, respectively.
subsequently, we incorporated irf <dig> mrna transcription by the active promoter, and to a lesser extent by the basal promoter state, as reactions  and , respectively.
reaction  considers the translation of irf <dig> mrna to produce irf <dig> protein.
irf <dig> protein is phosphorylated in reaction .
two phosphorylated irf <dig> proteins form a irf <dig> dimer in reaction .
finally, irf <dig> promoter inactivation, degradation of mrnas, isgf <dig>  ifn and irf <dig> proteins are represented by reactions  to , respectively:
<dig> ifn→fifn2ifn,
<dig> ifn→fisgf3isgf3+ifn,
<dig> ∅→fpapa,
<dig> pa→fmrnaamrna+pa,
<dig> ∅→fmrnabmrna,
<dig> mrna→firf7irf7+mrna,
<dig> irf7→firf7phosirf7phosp,
<dig> 2·irf7phosp→firf7dimerirf7dimer,
<dig> pa→fdpa∅,
<dig> ifn→fdifn∅,
<dig> mrna→fdmrna∅,
<dig> irf7dimer→fdirf7dimer∅,
<dig> isgf3→fdisgf3∅,
the reaction rates are given in table  <dig>
f
ifn
vifnifnnkifnn+ifnn
f
difn
k
difn·i
f
n
f
isgf3
k
isgf3·i
f
n
f
disgf3
k
disgf3·i
s
g
fpa
konisgf3·irf7dimerkai3·kai7+kai3·isgf3+kai7·irf7dimer+isgf3·irf7dimer
fdpa
k
off·p
a
fmrnaa
k
active·p
a
fmrnab
k
f
dmrna
k
dmrna·m
r
n
a
f
irf7
k
irf7·m
r
n
a
f
irf7phosp
k
dirf7·i
r
f
irf7dimer
k
irf7dimer·i
r
f7p
h
o
s
p·i
r
f7p
h
o
s
p
f
dirf7dimer
k
dirf7dimer·i
r
f7d
i
m
e
r
mathematical equations
to evaluate the deterministic precondition we consider the corresponding odes  to ):
<dig> difndt=fifn−fdifn
<dig> disgf3dt=fisgf3−fdisgf <dig>
<dig> dpadt=fpa−fdpa
<dig> dmrnadt=fmrnaa+fmrnab−fdmrna
<dig> dirf7dt=firf7−firf7phosp
<dig> dirf7phospdt=firf7phosp−2·firf7dimer
<dig> dirf7dimerdt=firf7dimer−fdirf7dimer
the model for the irf <dig> circuit can be obtained from biomodels database under reference model <dig>
fitting the stochastic irf <dig> gene expression model to experimental data
all parameters of the above described model for the irf <dig> gene circuit were fitted by using experimental data of irf protein expression.
in the model, a global quantity to describe the different forms of the irf <dig> protein was defined as follows:
<dig> irftotal=irf7+irf7phosp+irf7dimer,
and this variable was mapped to irf7†  flow cytometry measurements as follows:
<dig> irf7†=φirf7total,
where φ is a scaling factor.
using the experimental data, pdfs were build and the modes in the distributions were determined using the peakfinder function  <cit>  obtaining two elements of α =  <cit>  a.u.
, where t
n =  <dig> h. φ relates the values of fluorescence with the molecular count described by the mathematical model.
unfortunately, no calibration curve is provided with the data to calculate this parameter.
for this reason, multiple values were tested for φ obtaining consistent results, for illustrative purposes we report φ =  <dig>  we defined the deterministic precondition using eq.
, and assigning minimum and maximum acceptance ranges by setting βklow =  <dig>  and βkup =  <dig> , for k =  <dig>   <dig>  the allowed ranges for the parameter values are given in table  <dig>  the deterministic precondition was introduced in two different optimization strategies, random search, and genetic algorithms.
in both optimization strategies  <dig> realizations of the stochastic simulations were performed if the parameter set fulfilled the deterministic precondition.
v
ifn
n
k
ifn
k
difn
k
isgf3
k
disgf3
k
on
k
ai3
k
ai7
k
off
k
active
a
k
basal
b
k
dmrna
c
k
irf7
k
dirf7
k
irf7dimer
k
dirf7dimer
φ
a
k
basal was calculated to be at least one order of magnitude smaller than k
active
bdegradation rates for the mrna were calculated assuming a mrna half-life in the order of minutes  <cit>
cbased on the average translation rate in nih3t <dig> cells  <cit>
using the random search strategy, we tested  <dig> parameter sets from which less than the 1% passed the deterministic precondition and were stochastically evaluated.
the reduction of parameter values allowed us to efficiently find a set of parameter values that reproduced the experimental data .
a complete analysis of the performance of the use of the deterministic precondition with the random search algorithm is given in additional file  <dig>  to test the reproducibility of the estimated parameter values, the method was repeated  <dig> times obtaining well-defined parameter distributions that show the predominant values, see additional file  <dig>
then, we implemented the deterministic precondition using a genetic algorithm with adaptive population size.
here we implemented an initial population of  <dig> individuals for the first generation, and  <dig> subsequent generations with  <dig> individuals.
notice, that a large initial population of parameters was needed by the expected high rejection percentage of parameters by the deterministic precondition during the first generation.
as parameters for the algorithm we used ε= <dig>  and μ= <dig> .
our simulation results showed a constant decrease in the value of the objective function value during the generations, which indicates progress during fitting.
by the use of the deterministic precondition 99% of the parameters were rejected in the first generation and in the subsequent generations around 30% of the parameter values were rejected, improving in this way the efficiency of the algorithm .
the comparison between the experimental data and the model simulation distributions is given in fig.
<dig> showing a high degree of agreement.
in addition, to check the validity of the methodology further we fit the model to two additional flow cytometry measurements that describe the stimulation of the cell culture with 100u and 250u of ifn.
the respective results are given in additional file  <dig>  complete analysis of the performance of the use of the deterministic precondition in the genetic algorithm is given in additional file  <dig>
fig.
<dig> genetic algorithm performance.
the deterministic precondition was introduced in a genetic algorithm strategy.
the genetic algorithm was implemented using an adaptive population size using a population of  <dig> individuals for the first generation, and a population of  <dig> individuals for the subsequent generations.
as algorithm parameters we used μ =  <dig> , and ε =  <dig> .
in the plot a decrease in the objective function value during the generations in the ga is shown
fig.
<dig> comparison between the experimental data and stochastic simulations.
result of the fitting of the model with the experimental data using the genetic algorithm and the deterministic precondition.
in the plots, the y-axis represents the normalized cell count and the x-axis represents the fluorescence quantity  associated with the expression of the irf <dig> protein.
in gray we present the histograms that represent the experimental data, in red the pdf from the stochastic simulations
irf <dig> temporal dynamics
the simple model of irf <dig> gene expression described above is sufficient to explain irf <dig> bimodality.
using the optimized parameter values given in table  <dig> and the initial conditions given in table  <dig>  the stochastic temporal dynamics of the model were simulated and are given in fig.
<dig>  in fig.
7-a it can be seen that the ifn concentration evolves to a steady state.
subsequently, the first affected variable is isgf <dig> that equally evolves towards a stable steady state, see fig.
7-b.
figure 7-c shows the irf <dig> promoter dynamics exhibiting transitions between active/basal states.
this promoter state transition is a characteristic of genes with regulated expression  <cit> .
subsequently, irf <dig> mrna expression displays a pattern that is affected by this stochastic switching between two possible states, one with basal expression and the other with active expression, see fig.
7-d. for single-cell trajectories, irf <dig> protein expression shows a switch-like expression.
for the whole population of those trajectories bimodality is observed , the same stands for the different forms of the irf <dig> protein, phosphorylated  dimer , and the total amount of irf <dig> proteins .
the system’s steady states are given in table  <dig>
fig.
<dig> time courses in the irf <dig> circuit after ifn stimulation.
the temporal dynamics of the promoter, mrna and irf <dig> protein dynamics were obtained after stochastic simulations.
a representative trajectory that presents a single cell dynamics is given by the red lines.
the stochastic simulations were repeated  <dig> times using the same initial condition obtaining the histograms that represent the cell population.
a time courses of ifn showing the evolution to a steady state, the same is observed in b for the isgf <dig>  c irf <dig> promoter shows the transition between the two possible states off/on.
d irf <dig> mrna showing a basal expression state and a state of active expression.
bimodality is observed in e for the irf <dig> protein, the same bimodal behavior was observed in f for the irf <dig> phosphorylated protein, in g for the irf <dig> protein dimer, and in h for all forms of the irf <dig> protein
ifn
i
s
g
pa
mrna
i
r
i
r
f7p
h
o
s
p
i
r
f7d
i
m
e
r
ifn
i
s
g
pa
mrna
i
r
i
r
f7p
h
o
s
p
i
r
f7d
i
m
e
r
i
r
f7t
o
t
a
l
discussion
the promise of systems biology is to achieve a quantitative understanding of the molecular processes in the cell with the aid of computational models.
however, a bottleneck is the availability of reliable parameter values needed in those models.
often, it is very difficult or even impossible to measure all of these parameters.
for deterministic models, this problem has been well tackled by the development of efficient methods for parameter estimation.
contrarily, for stochastic dynamics the landscape of methods for parameter estimation is poorer and still under development.
for this reason, innovative and efficient algorithms are needed to fit and validate realistic stochastic models.
during the last years important advances have been achieved in the field of parameter inference for stochastic models.
on one hand, strategies that involve mathematical procedures of moment estimation have been suggested .
on the other, bayesian methods that test multiple parameters values to find approximated solutions that represent the data have been successfully implemented  <cit> .
currently, the main problems observed in most methods include the computational cost, the accuracy of the obtained solution, and its potential to be implemented in large and non-linear systems.
different strategies have been suggested to improve the accuracy and alleviate the computational performance  <cit> .
in our method, we tackled the computational cost by introducing a deterministic precondition that works as a by-pass in the algorithm avoiding large amounts of unnecessary stochastic simulations.
the concept of reducing the parameter space by introducing a precondition defined by the experimental data has been suggested previously by hori et al.
<cit> .
in their method, a small order linear model is used to optimize the experimental data by finding the root of a lyapunov equation.
in contrast, we use a deterministic precondition.
both approaches have advantages and disadvantages.
thus, hori’s method is constrained by the need to find the lyapunov equation.
our method is based on steady state calculation and monte carlo simulations that are standard and well-known methodologies used in systems biology.
recently a new algorithm developed by lillacci et al., has been shown to significantly reduce the computational cost and achieve a high accuracy fit for stochastic models and flow cytometry data  <cit> .
this method uses the kolmogorov distance as a metric to calculate the difference between model simulation and experimental data.
by choosing this metric, it is possible to estimate the minimal number of simulations needed to compare experimental and model simulations under a certain tolerance value.
this estimated value decreases as the number of experimental data increases.
in typical flow cytometry experiments the number of measured cells is in the order of tens of thousands and this number in lillacci’s algorithm is translated in a reduction of at least one order of magnitude in the number of required stochastic simulations.
this algorithm has been applied to a model with  <dig> reactions and  <dig> free parameters obtaining a good fit to flow cytometry experimental data that reproduces bimodal dynamics.
comparing our method with lillacci’s algorithm it is important to point out that both methods tackle the computational cost issue in two fundamentally different ways.
lillacci’s algorithm minimizes the number of needed stochastic simulations, whereas our algorithm minimizes the number of parameters sets evaluated with stochastic dynamics.
a powerful new algorithm may be the result of combining both methodologies.
it is well known from non-linear dynamics that the model architecture and parameter values determine the behavior of the system.
hence, more complex stable dynamics such as limit cycles, and higher-multistability  can be obtained.
unstable systems are usually not of interest in the biological context.
our method was designed and tested for monostable and bistable systems, all other cases being rejected by the deterministic precondition.
the case of limit cycles can be better approached by existing parameter estimation methods for single time-course data .
the case of systems with higher-multistability can in principle be treated by our method, but the application is limited by the costly need to find multiple stable steady-states in the system.
finally, it is relevant to mention that comparing stable states of the ode model with modes of the measured flow cytometry distribution in some models can present some inherent problems.
bimodal fluorescence distributions may have no relation to bistability, e.g.
when the systems has large differences in transitions rates for the promoter states  <cit> .
our method was implemented using two well-established optimization strategies, random search and genetic algorithms.
using a random search algorithm we observed that less than 1% of the population of parameters are subject to stochastic evaluations, this reduction in the tested parameters allowed us to explore a larger proportion of the parameter space, which is especially relevant for systems with multiple unknown parameters and no initial parameter guesses.
on the other hand, using genetic algorithms we achieved a convergence to a minimal value in the optimization function after a few generations .
additionally, during each generation in the genetic algorithm a large percentage of the total evaluated parameters is rejected by the deterministic precondition.
this percentage is dependent on the parameter of the algorithm .
in both strategies the introduction of the deterministic precondition significantly improved the parameter estimation process.
moreover, the implementation of deterministic precondition is not restricted to random search and genetic algorithms, it can also be implemented in other optimization algorithms, e.g.
particle swarm or simulated annealing.
a complete analysis of the performance of using the deterministic condition is given in additional file  <dig>
the accuracy of the method is given by the agreement between experimental data and simulations.
as can be observed in fig.
<dig> a good fit was found, albeit not a perfect one.
we observed that even increasing the number of tested parameter values during the random search and/or increasing the number of generations during the genetic algorithm did not improve the fit.
for this reason, we consider that the differences between the model and the experimental data might be explained by the fact that our system is a highly simplified model that was built by lumping some steps in the biological system.
our team is working to fit more complex models in a future publication.
to test the reproducibility of the obtained parameter values, the method was repeated  <dig> times and distributions of the obtained parameters were built.
as can be observed in additional file  <dig>  when a parameter is identifiable this is reflected in a narrow distribution, on the other hand, when a parameter is non-identifiable this is reflected in a wide distribution.
in the example given by reactions  to  we observed that most parameters are not identifiable.
thus, non-identifiability is expected in the parameters contained in reactions  and .
those reactions have opposite effects on the ifn dynamics.
for this reason, multiple combinations of values in parameters v
ifn, k
ifn and k
difn have similar effects on the overall system dynamics.
contrarily, the parameters involved in the production of the mrna  are better confined.
rand et al.
described bimodal gene expression of irf <dig> after ifn stimulation.
this phenomenon was explained at the cell population level by the effects of the ifn paracrine response  <cit> .
here we observed that bimodality also can be explained in absence of paracrine response and only taking into account the molecular mechanism responsible for the irf <dig> promoter dynamics.
bistability in gene-expression circuits is commonly associated with the switching between active/basal states in the gene promoter.
in most cases, dna cooperativity in the promoter is the basis of promoter-state switching  <cit> .
however, for type-i ifn responses a promoter activation in a cooperative manner has recently been discarded  <cit> .
taking the recent literature into account, we developed a model on the basis of a positive feedback loop circuit comprising the independent activation of isre and irfe elements by one isgf <dig> molecule and one irf <dig> dimer, respectively .
this model has the needed and sufficient elements to sustain a bistable system   <cit> .
additionally, multiple model structures testing different biological scenarios including: additional feedback loops, additional intermediate elements, and promoter cooperative activation were tested, obtaining that the presented model reproduces the experimental data best.
our simulation results agree with the experimental data obtained from flow cytometry  <cit> .
additional characteristics observed in the experimental data are also reproduced by the model, such as the basal expression of irf <dig> without ifn stimulation  <cit> , see additional file  <dig>  the number of molecules in the active state of the system  agrees with the order of magnitude reported for mammal mrnas  and the range of protein concentration   <cit> .
CONCLUSIONS
here, we present a method to fit stochastic models to experimental data.
the method is based on the comparison of distributions.
the central idea of the method is to use a deterministic precondition that is defined by the experimental data as a filter avoiding large amounts of costly stochastic simulations.
using this idea, the number of parameters evaluated under stochastic dynamics is reduced, resulting in a significant improvement in the performance of the algorithm.
as a case study, we used a model of irf <dig> gene expression investigating the origin of bimodality in its dynamics upon ifn stimulation.
our results allowed us to conclude that a circuit with irf <dig> promoter activation by one irf <dig> dimer and one isgf <dig> molecule is sufficient to explain the observed bimodal dynamics.
additional files
additional file  <dig> pseudocode for the genetic algorithm with the deterministic precondition.
additional file with the pseudocode for the genetic algorithm.
additional file  <dig> identifying parameters for a constitutive gene expression circuit with in-silico data.
additional file with an example where the method is applied to a simulation case with synthetic data.
the file contains figures a1-a <dig>  and tables a <dig> and a <dig>
additional file  <dig> random search strategy.
additional file with figure a <dig>  results obtained by using the random search strategy.
additional file  <dig> algorithm’s speed and performance analysis.
additional file with an analysis of the algorithm’s speed and performance.
the file contains figure a <dig>  and tables a <dig> and a <dig>
additional file  <dig> testing the reproducibility of the selected parameter values.
additional file containing figure a <dig> that shows the parameter distributions.
additional file  <dig> fitting the model to different doses of ifn.
additional file with figure a <dig> and a <dig> that shows the fits of the model to flow cytometry data describing the stimulation of cells with 100u and 250u of ifn.
additional file  <dig> model dynamics without ifn stimulation.
additional file with figure a <dig> that shows the time courses for the model in absence of ifn stimulation.
abbreviations
auarbitrary units
bacbacterial artificial chromosome
dpdeterministic precondition
gagenetic algorithm
ifninterferon
isginterferon stimulated genes
odeordinary differential equation
pdfprobability density function
rsrandom search
sssstable steady state

 computer models contribute to the integrative understanding of complex molecular processes in the cell. the most commonly used approach is deterministic modeling based on ordinary differential equations. when the studied system comprises species with a low molecular count, stochastic formalisms, e.g. gillespie’s algorithm that simulates trajectories and uses discrete molecule numbers are more appropriate. in recent years, the use of stochastic models has substantially increased. considering stochasticity in biological systems has changed the quantitative and qualitative understanding obtained by previous deterministic models. examples of biological phenomena discovered by stochastic modeling include gene expression in burst-like patterns, productive or latent cell decision after hiv-infection, and the presence of oscillatory behavior induced by noise. however, in contrast to the plethora of methods used to fit and analyze deterministic models, there are only very limited sets of methods available to do the same with stochastic models. thus, only recently, methods for parameter fitting of stochastic models have been developed and are still under development. parameter estimation methods for stochastic models have been so far designed for single time course data. for experiments leading to data distributions, they include the moment closure, and the comparison of experimental and simulation distributions. despite the current efforts, there are still major problems that limit the full applicability of those methods. the main drawback is the high computational execution time. to solve this problem novel efficient strategies to fit stochastic models are needed. in the following we introduce a new approach that is a variant of existing methods based on the comparison of distributions. differently to the existing methods we use the experimental data to define a condition that the model must fulfill in the deterministic regime. this condition tests, if the deterministic steady states are in the close proximity to the modes of the experimental distribution. only if the model fulfills this condition, the parameter set is evaluated in the stochastic regime. in this way, the algorithm directs the evaluation of parameter sets towards regions in parameter space with high potential to reproduce the experimental data. additionally, this method applies to non-linear models and complex experimental data distributions. to test the potential of this new method, we selected an open scientific question and real experimental data. we investigated the molecular mechanisms responsible for the experimentally observed bimodal dynamics of irf expression after interferon stimulation. it is well documented that in isogenic cell populations not all cells produce an antiviral response when infected by a virus. the mechanisms behind this heterogeneity have been related to stochastic events in the signaling pathways responsible to elicit the antiviral response. ifn is a cytokine that activates the jak-stat signaling pathway in virus-infected cells. the jak-stat signaling pathway induces the translocation of isgf into the nucleus to directly activate the transcription of a set of several hundred ifn-stimulated genes. irf is an isg with a central role in the immune response. recently, it was observed that after ifn stimulation murine fibroblasts show a switch-like pattern of irf expression, which is reflected at the population level as bimodality. we developed a model to describe the observed bimodal dynamics of irf expression after ifn stimulation. we hypothesize that the binding of irf and isgf to the irf promoter is the mechanism responsible for this bimodal behavior. our simulation results reproduced irf switch-like dynamics at single cell level, and bimodality was achieved at the population level. hereby, we used the newly developed method to fit the model and correctly reproduce the experimental data. our results allow us to conclude that the irf promoter activation by the combination of irf dimers and isgf is sufficient to quantitatively explain the observed irf bimodal dynamics. methods experimental data published experimental data were produced by rand et al.. the data described the expression of irf after ifn- β stimulation in a population of murine nih3t fibroblasts. rand’s experiments were done in the following way: first, cells were transfected with a bac containing irf and reporter mcherry genes fused, subsequently cultures were treated with different concentrations of murine ifn- β. for illustrative purposes we selected the treatment with 150u of ifn- β, a concentration where bimodality is prevalent. then, fluorescence in the cultures was monitored using flow cytometry at different time points during hours. numerical methods in the deterministic regime the model was simulated using symbolic methods in matlab and/or using the lsoda algorithm in copasi. stability was calculated by making the right-hand side of the differential equations equal to zero and subsequently by finding the roots of the system using function solve in matlab. eigenvalues were calculated using function eig in matlab. for solving the model under stochastic dynamics we used the gibson and bruck algorithm coded in copasi. the random search and genetic algorithm were coded in matlab. raw experimental data was analyzed using the function fcs data reader coded in matlab. the modes in the distributions were calculated using the function peakfinder coded in matlab. the source code of the project can be accessed via: https://sourceforge.net/projects/irf7-bimodaldynamics/.  new algorithm to fit stochastic biological models comparing experimental and simulation distributions the measurements of fluorescence were made comparable with the corresponding observable chemical species s o in the model by a function h that maps state s o to observation s † as follows: s†=h), at each measurement time point t i, i= …,n. subsequently, a process to compare distributions was developed based on. first, considering a specific set of parameter values θ={θ …,θ d}, we performed ns repetitions of the stochastic simulations s={s …,s ns}. the total of those stochastic simulations were used to build histograms with a fixed number of bins l for each t i. subsequently, the probability density function of the simulations p s,θ,b l) were computed using the center of each bin b l, for l= …,l, and the normalized form of the histograms. on the other hand, having nm repetitions of single cell experimental data m={m …,m nm} histograms and pdfs p e,b l) were built. even though, complex formulas can be use to calculate the number of bins in an histogram, here l was kept constant for p s and p e and was calculated as l=nm, this is a simple approach used by default by most programming languages. finally, an approach using the difference of squares was used as a metric to calculate the distance between p s,θ,b l) and p e,b l), as follows: f=∑i=1n∑l=1l,bl)−ps,θ,bl)) algorithm implementation the objective function f can be used to calculate a parameter estimate θ^=argminθf. as optimization usually requires plenty of function evaluations and each evaluation of f requires ns stochastic simulations, a direct optimization strategy is computationally unfeasible even for simple models. therefore, we introduce a new strategy that selects good candidate parameters and only performs the stochastic simulations for these parameters. in large systems where fluctuations can be discarded, the stochastic system can be reduced to the deterministic one. for this reason, in most cases the deterministic dynamics can be associated with a measure of central tendency in the pdfs obtained after the stochastic simulations. using this reasoning we will introduce the strategy to efficiently estimate parameters for stochastic models making use of deterministic dynamics as an initial indicator. assuming that the experimental data is in equilibrium at the last measurement point t n, we used the modes from p e,b l) as a central measure of tendency, denoting the modes by α=,…,α q), with q being the total number of modes. then, using the ode version of the model ) we determine its stability. here we tested whether the system has stable steady states and denote them with x∗=x1∗,…,xss∗, being s s the number of stable steady states. if the system has no stable steady state, it holds s s= and the vector x ∗ is empty. the calculation of steady states can either be performed by solving y˙= where y˙ is the right hand side of the ode system with x replaced by y. solving this equation can be impossible for large systems. alternatively, one can numerically solve the ode systems until the flux is zero for all components. varying the initial conditions leads to the different steady states. the stability of the system was calculated after determining the sign of the eigenvalues, for stable steady states. having x ∗ and α we introduce the deterministic precondition: βklowxk∗≤αk≤βkupxk∗,for1≤k≤q,and0≤βklow≤1and1≤βkup. equation means that the number of modes equals the number of stable steady states: s s=q, and that each of the stable steady states is close to its corresponding mode. only, if this is fulfilled, we calculate stochastic simulations and evaluate the objective function f. if the precondition is not fulfilled, we do not carry out stochastic simulations. in this case, we assign a high objective function value to direct the parameter search towards parameters that pass the deterministic precondition. a graphical representation of the deterministic precondition is given in fig. fig. the deterministic precondition. left side the deterministic dynamics of the system can be seen, here the observable variable evolves to a single stable steady state x1∗. this steady state can be calculated by finding the point in the system where the net flux is equal to zero x˙= at the right hand side of the plot, the pdf of the experimental data is shown. this pdf is rotated and the x-axis shows the count and the y-axis shows the bin values with units of molecules/cell. from both graphs it can be seen that the deterministic precondition is valid when x1∗ is inside the range defined by β1lowx1∗ and β1upx1∗ at this point, it is important to notice that the deterministic precondition is not stating that the stable steady states form the ode system are equal to the mean obtained by solving the stochastic system. rather, the deterministic precondition tests if the deterministic steady state lies around a certain range. fcond=f,βklowxk∗≤αk≤βkupxk∗−−′′condition passed′′∞,else –′′condition not passed′′ our extended objective function with the deterministic check reads as: and the estimate is defined as: θ^cond=argminθfcond. as previously discussed, our method is based on the central assumption that the modes in the experimental data are related to the stable steady states in the deterministic mathematical model. if this assumption should fail, either of the following could happen: a) the method accepts a parameter that passed the deterministic precondition but does not show a good agreement of the distributions. if it does not show good agreement, the objective function f will show a high value which means that this parameter does not lead to a good fit. therefore, this case leads to a loss of computational time, but not to wrong results and it overall is still less expensive than evaluating every parameter set. b) we reject a parameter that fails the deterministic precondition but would have lead to a good fit. in this case, we lose indeed a good parameter. this shows that our method is obviously a heuristic strategy that does not guarantee to find the exact global minimum, but this is true for almost any optimization strategy. the algorithmic steps are depicted in fig. fig. algorithm to fit stochastic models to experimental data. the algorithm solves the model in the deterministic and stochastic regimes. a condition observed in the experimental data is defined. a set of parameter values is evaluated in the deterministic regime to test if the model reproduces this condition. if the condition is met the stochastic simulation is performed. otherwise, the parameter values are rejected. the pdf obtained from the experimental data is compared with the pdf obtained after running the stochastic simulations. this comparison is made using a difference of squares as an objective function. this process is repeated until evaluating a total number of parameter sets or after a termination criterion is met. the parameter set that best reproduces the experimental data is given by the minimum value obtained after the iterative evaluation of the objective function parameter estimation methods the deterministic precondition was first implemented using a random search algorithm. random search is a global optimization strategy that tests random combinations of parameter values. the successful output in this method is dependent on the total number of evaluated parameter sets. the more evaluated parameter values, the higher the probability to find the global minimum. the pseudo-code for the random search is given in algorithm the second implementation of the deterministic precondition was using the more directed genetic algorithm. genetic algorithm mimics evolution and is based on reproduction and selection. this algorithm is made of a population of individuals, and each contains a genome that is defined by the number of parameters to optimize. the individuals are ranked after solving the objective function, and a population of parental individuals is selected according to an elitism rate. new individuals are generated by pairing and recombining the parental genomes. variability is introduced in the population by adding mutations in the new individuals according to a given mutation rate. by the continuous process of selecting the best parameters after each generation, the algorithm evolves towards the regions in the parameter space that reduce the values in the objective function. the pseudo-code for the genetic algorithm is given in additional file identifying parameters for a constitutive gene expression circuit with in-silico data to prove the functionality and benefits of the new proposed algorithm, we first applied it to a simple example where the parameters are known a priori. the model describes the stochastic dynamics of two variables, protein and mrna of a gene with constitutive expression, a graphical representation of the constitutive gene expression circuit is given in fig. 3-a. the model is described by the following reactions: fig. parameter estimation using a constitutive gene expression circuit and in-silico data. a diagram of a gene with constitutive expression. the mrna is produced at constant rate θ and is degraded at constant rate θ the proteins are produced at constant rate θ and are degraded with θ b comparison between in-silico data and model dynamics using the true parameters θ o= and the estimated parameters θ^=. a distribution with the in-silico data is given in grey. in red is given the model dynamics. c comparison between a priori and a posteriori parameter distributions. the true parameters are represented by the intersection of the red lines. in the priori distribution it can be observed that parameters are randomly distributed in the parameter space. the posteriori distributions were calculated by running independent random searches with parameters each, and by plotting the best parameter value selected by the algorithm. in the posteriori distribution it is plotted only the parameters that passed the deterministic precondition, and it can be observed that those parameter estimates are in a region close to the true parameter values ∅→θ1mrna, mrna→θ2∅, mrna→θ3protein, protein→θ4∅. to estimate the parameter values for this model we carried out the following procedure: first, we generated the in-silico data by selecting the following true parameter values θ o= and running independent stochastic simulations from which pdfs were built. four observable time points were defined at min, min, min and min. from the pdfs a mode was obtained at α = a.u., where t n= min. then, we defined the deterministic precondition using eq., and assigning minimum and maximum acceptance ranges by setting β1low = and β1up =, respectively. a deeper analysis of the stability of the constitutive gene expression circuit is given in additional file to build simulated pdfs we used repetitions of the stochastic model, this number was empirically calculated according to additional file for illustrative purposes, we assumed that the values for the parameters responsible for the mrna transcription and protein translation were unknown, and the new algorithm was used to estimate those parameters. the deterministic precondition was applied using the rs strategy obtaining that the algorithm only evaluates stochastic dynamics in % of the tested parameter values, reducing in this way the total simulation time. additionally, the complete algorithm was repeated times and histograms of the estimated parameter were computed to determine whether they are close to θ. as can be observed in fig. 3-c the deterministic precondition reduces the evaluation of different parameters under stochastic dynamics by selecting only those parameters that are in a well-defined region in the proximity of θ. for this model, the main benefits of using the deterministic precondition was the reduction of the number of parameters evaluated under stochastic dynamics. this rejection of parameters was made in an area outside the true parameter values, and hence no difference in accuracy is expected in comparison with a method without using the deterministic precondition. a complete description of the analysis of the performance, accuracy and error for this example is given in additional file the model for the constitutive gene expression circuit can be obtained from biomodels database under reference model mathematical model for irf expression dynamics subsequently we applied our algorithm to a real problem with flow cytometry data. here we studied the dynamics of murine irf gene expression upon ifn stimulation. for this reason, we developed a model that comprises known key components and feedback mechanisms. the overall system describes the active irf promoter by the binding of irf dimer and isgf to the dna binding sites isre and irfe, respectively, the transcription and translation of irf and its subsequent phosphorylation and dimerization. irf protein binding to the irfe binding site in the promoter results in the production of more irf protein, constituting a positive feedback loop. a graphical representation of the irf gene expression dynamics is given in fig. fig. irf gene expression circuit. a in viral infected cells irf gene expression is induced after ifn stimulation by the jak-stat signaling pathway. the irf promoter activation is governed by the binding of irf dimer and isgf to the promoter dna binding sites. irf promoter activation leads to the transcription of irf mrna and subsequent its translation to produce the irf protein. notice that the irf production and subsequently phosphorylation and dimerization leads to the binding to its own dna binding site, resulting in the production of more irf making in this way a positive feedback loop. b a simplified system representing the irf gene expression dynamics is given. here, the promoter transitions between active/basal states, the gene transcription and translation processes are represented as solid black lines as well as the feedback loop in the system chemical reactions the developed model consists of species and reactions to ). in the model, reaction describes the ifn activation of the jak-stat signaling pathway. for the description of this reaction and according to a saturable function was used. downstream reactions of the pathway were lumped in reaction, so that the same variable was used to describe the output of the jak-stat signaling pathway, namely isgf in this highly simplified model. reaction describes the irf promoter activation by the binding of irf dimer and isgf to isre and irfe, respectively. subsequently, we incorporated irf mrna transcription by the active promoter, and to a lesser extent by the basal promoter state, as reactions and, respectively. reaction considers the translation of irf mrna to produce irf protein. irf protein is phosphorylated in reaction. two phosphorylated irf proteins form a irf dimer in reaction. finally, irf promoter inactivation, degradation of mrnas, isgf ifn and irf proteins are represented by reactions to, respectively: ifn→fifn2ifn, ifn→fisgf3isgf3+ifn, ∅→fpapa, pa→fmrnaamrna+pa, ∅→fmrnabmrna, mrna→firf7irf7+mrna, irf7→firf7phosirf7phosp, 2·irf7phosp→firf7dimerirf7dimer, pa→fdpa∅, ifn→fdifn∅, mrna→fdmrna∅, irf7dimer→fdirf7dimer∅, isgf3→fdisgf3∅, the reaction rates are given in table f ifn vifnifnnkifnn+ifnn f difn k difn·i f n f isgf3 k isgf3·i f n f disgf3 k disgf3·i s g fpa konisgf3·irf7dimerkai3·kai7+kai3·isgf3+kai7·irf7dimer+isgf3·irf7dimer fdpa k off·p a fmrnaa k active·p a fmrnab k f dmrna k dmrna·m r n a f irf7 k irf7·m r n a f irf7phosp k dirf7·i r f irf7dimer k irf7dimer·i r f7p h o s p·i r f7p h o s p f dirf7dimer k dirf7dimer·i r f7d i m e r mathematical equations to evaluate the deterministic precondition we consider the corresponding odes to ): difndt=fifn−fdifn disgf3dt=fisgf3−fdisgf dpadt=fpa−fdpa dmrnadt=fmrnaa+fmrnab−fdmrna dirf7dt=firf7−firf7phosp dirf7phospdt=firf7phosp−2·firf7dimer dirf7dimerdt=firf7dimer−fdirf7dimer the model for the irf circuit can be obtained from biomodels database under reference model fitting the stochastic irf gene expression model to experimental data all parameters of the above described model for the irf gene circuit were fitted by using experimental data of irf protein expression. in the model, a global quantity to describe the different forms of the irf protein was defined as follows: irftotal=irf7+irf7phosp+irf7dimer, and this variable was mapped to irf7† flow cytometry measurements as follows: irf7†=φirf7total, where φ is a scaling factor. using the experimental data, pdfs were build and the modes in the distributions were determined using the peakfinder function obtaining two elements of α = a.u., where t n = h. φ relates the values of fluorescence with the molecular count described by the mathematical model. unfortunately, no calibration curve is provided with the data to calculate this parameter. for this reason, multiple values were tested for φ obtaining consistent results, for illustrative purposes we report φ = we defined the deterministic precondition using eq., and assigning minimum and maximum acceptance ranges by setting βklow = and βkup =, for k = the allowed ranges for the parameter values are given in table the deterministic precondition was introduced in two different optimization strategies, random search, and genetic algorithms. in both optimization strategies realizations of the stochastic simulations were performed if the parameter set fulfilled the deterministic precondition. v ifn n k ifn k difn k isgf3 k disgf3 k on k ai3 k ai7 k off k active a k basal b k dmrna c k irf7 k dirf7 k irf7dimer k dirf7dimer φ a k basal was calculated to be at least one order of magnitude smaller than k active bdegradation rates for the mrna were calculated assuming a mrna half-life in the order of minutes cbased on the average translation rate in nih3t cells using the random search strategy, we tested parameter sets from which less than the 1% passed the deterministic precondition and were stochastically evaluated. the reduction of parameter values allowed us to efficiently find a set of parameter values that reproduced the experimental data. a complete analysis of the performance of the use of the deterministic precondition with the random search algorithm is given in additional file to test the reproducibility of the estimated parameter values, the method was repeated times obtaining well-defined parameter distributions that show the predominant values, see additional file then, we implemented the deterministic precondition using a genetic algorithm with adaptive population size. here we implemented an initial population of individuals for the first generation, and subsequent generations with individuals. notice, that a large initial population of parameters was needed by the expected high rejection percentage of parameters by the deterministic precondition during the first generation. as parameters for the algorithm we used ε= and μ=. our simulation results showed a constant decrease in the value of the objective function value during the generations, which indicates progress during fitting. by the use of the deterministic precondition 99% of the parameters were rejected in the first generation and in the subsequent generations around 30% of the parameter values were rejected, improving in this way the efficiency of the algorithm. the comparison between the experimental data and the model simulation distributions is given in fig. showing a high degree of agreement. in addition, to check the validity of the methodology further we fit the model to two additional flow cytometry measurements that describe the stimulation of the cell culture with 100u and 250u of ifn. the respective results are given in additional file complete analysis of the performance of the use of the deterministic precondition in the genetic algorithm is given in additional file fig. genetic algorithm performance. the deterministic precondition was introduced in a genetic algorithm strategy. the genetic algorithm was implemented using an adaptive population size using a population of individuals for the first generation, and a population of individuals for the subsequent generations. as algorithm parameters we used μ =, and ε =. in the plot a decrease in the objective function value during the generations in the ga is shown fig. comparison between the experimental data and stochastic simulations. result of the fitting of the model with the experimental data using the genetic algorithm and the deterministic precondition. in the plots, the y-axis represents the normalized cell count and the x-axis represents the fluorescence quantity associated with the expression of the irf protein. in gray we present the histograms that represent the experimental data, in red the pdf from the stochastic simulations irf temporal dynamics the simple model of irf gene expression described above is sufficient to explain irf bimodality. using the optimized parameter values given in table and the initial conditions given in table the stochastic temporal dynamics of the model were simulated and are given in fig. in fig. 7-a it can be seen that the ifn concentration evolves to a steady state. subsequently, the first affected variable is isgf that equally evolves towards a stable steady state, see fig. 7-b. figure 7-c shows the irf promoter dynamics exhibiting transitions between active/basal states. this promoter state transition is a characteristic of genes with regulated expression. subsequently, irf mrna expression displays a pattern that is affected by this stochastic switching between two possible states, one with basal expression and the other with active expression, see fig. 7-d. for single-cell trajectories, irf protein expression shows a switch-like expression. for the whole population of those trajectories bimodality is observed, the same stands for the different forms of the irf protein, phosphorylated dimer, and the total amount of irf proteins. the system’s steady states are given in table fig. time courses in the irf circuit after ifn stimulation. the temporal dynamics of the promoter, mrna and irf protein dynamics were obtained after stochastic simulations. a representative trajectory that presents a single cell dynamics is given by the red lines. the stochastic simulations were repeated times using the same initial condition obtaining the histograms that represent the cell population. a time courses of ifn showing the evolution to a steady state, the same is observed in b for the isgf c irf promoter shows the transition between the two possible states off/on. d irf mrna showing a basal expression state and a state of active expression. bimodality is observed in e for the irf protein, the same bimodal behavior was observed in f for the irf phosphorylated protein, in g for the irf protein dimer, and in h for all forms of the irf protein ifn i s g pa mrna i r i r f7p h o s p i r f7d i m e r ifn i s g pa mrna i r i r f7p h o s p i r f7d i m e r i r f7t o t a l discussion the promise of systems biology is to achieve a quantitative understanding of the molecular processes in the cell with the aid of computational models. however, a bottleneck is the availability of reliable parameter values needed in those models. often, it is very difficult or even impossible to measure all of these parameters. for deterministic models, this problem has been well tackled by the development of efficient methods for parameter estimation. contrarily, for stochastic dynamics the landscape of methods for parameter estimation is poorer and still under development. for this reason, innovative and efficient algorithms are needed to fit and validate realistic stochastic models. during the last years important advances have been achieved in the field of parameter inference for stochastic models. on one hand, strategies that involve mathematical procedures of moment estimation have been suggested. on the other, bayesian methods that test multiple parameters values to find approximated solutions that represent the data have been successfully implemented. currently, the main problems observed in most methods include the computational cost, the accuracy of the obtained solution, and its potential to be implemented in large and non-linear systems. different strategies have been suggested to improve the accuracy and alleviate the computational performance. in our method, we tackled the computational cost by introducing a deterministic precondition that works as a by-pass in the algorithm avoiding large amounts of unnecessary stochastic simulations. the concept of reducing the parameter space by introducing a precondition defined by the experimental data has been suggested previously by hori et al.. in their method, a small order linear model is used to optimize the experimental data by finding the root of a lyapunov equation. in contrast, we use a deterministic precondition. both approaches have advantages and disadvantages. thus, hori’s method is constrained by the need to find the lyapunov equation. our method is based on steady state calculation and monte carlo simulations that are standard and well-known methodologies used in systems biology. recently a new algorithm developed by lillacci et al., has been shown to significantly reduce the computational cost and achieve a high accuracy fit for stochastic models and flow cytometry data. this method uses the kolmogorov distance as a metric to calculate the difference between model simulation and experimental data. by choosing this metric, it is possible to estimate the minimal number of simulations needed to compare experimental and model simulations under a certain tolerance value. this estimated value decreases as the number of experimental data increases. in typical flow cytometry experiments the number of measured cells is in the order of tens of thousands and this number in lillacci’s algorithm is translated in a reduction of at least one order of magnitude in the number of required stochastic simulations. this algorithm has been applied to a model with reactions and free parameters obtaining a good fit to flow cytometry experimental data that reproduces bimodal dynamics. comparing our method with lillacci’s algorithm it is important to point out that both methods tackle the computational cost issue in two fundamentally different ways. lillacci’s algorithm minimizes the number of needed stochastic simulations, whereas our algorithm minimizes the number of parameters sets evaluated with stochastic dynamics. a powerful new algorithm may be the result of combining both methodologies. it is well known from non-linear dynamics that the model architecture and parameter values determine the behavior of the system. hence, more complex stable dynamics such as limit cycles, and higher-multistability can be obtained. unstable systems are usually not of interest in the biological context. our method was designed and tested for monostable and bistable systems, all other cases being rejected by the deterministic precondition. the case of limit cycles can be better approached by existing parameter estimation methods for single time-course data. the case of systems with higher-multistability can in principle be treated by our method, but the application is limited by the costly need to find multiple stable steady-states in the system. finally, it is relevant to mention that comparing stable states of the ode model with modes of the measured flow cytometry distribution in some models can present some inherent problems. bimodal fluorescence distributions may have no relation to bistability, e.g. when the systems has large differences in transitions rates for the promoter states. our method was implemented using two well-established optimization strategies, random search and genetic algorithms. using a random search algorithm we observed that less than 1% of the population of parameters are subject to stochastic evaluations, this reduction in the tested parameters allowed us to explore a larger proportion of the parameter space, which is especially relevant for systems with multiple unknown parameters and no initial parameter guesses. on the other hand, using genetic algorithms we achieved a convergence to a minimal value in the optimization function after a few generations. additionally, during each generation in the genetic algorithm a large percentage of the total evaluated parameters is rejected by the deterministic precondition. this percentage is dependent on the parameter of the algorithm. in both strategies the introduction of the deterministic precondition significantly improved the parameter estimation process. moreover, the implementation of deterministic precondition is not restricted to random search and genetic algorithms, it can also be implemented in other optimization algorithms, e.g. particle swarm or simulated annealing. a complete analysis of the performance of using the deterministic condition is given in additional file the accuracy of the method is given by the agreement between experimental data and simulations. as can be observed in fig. a good fit was found, albeit not a perfect one. we observed that even increasing the number of tested parameter values during the random search and/or increasing the number of generations during the genetic algorithm did not improve the fit. for this reason, we consider that the differences between the model and the experimental data might be explained by the fact that our system is a highly simplified model that was built by lumping some steps in the biological system. our team is working to fit more complex models in a future publication. to test the reproducibility of the obtained parameter values, the method was repeated times and distributions of the obtained parameters were built. as can be observed in additional file when a parameter is identifiable this is reflected in a narrow distribution, on the other hand, when a parameter is non-identifiable this is reflected in a wide distribution. in the example given by reactions to we observed that most parameters are not identifiable. thus, non-identifiability is expected in the parameters contained in reactions and. those reactions have opposite effects on the ifn dynamics. for this reason, multiple combinations of values in parameters v ifn, k ifn and k difn have similar effects on the overall system dynamics. contrarily, the parameters involved in the production of the mrna are better confined. rand et al. described bimodal gene expression of irf after ifn stimulation. this phenomenon was explained at the cell population level by the effects of the ifn paracrine response. here we observed that bimodality also can be explained in absence of paracrine response and only taking into account the molecular mechanism responsible for the irf promoter dynamics. bistability in gene-expression circuits is commonly associated with the switching between active/basal states in the gene promoter. in most cases, dna cooperativity in the promoter is the basis of promoter-state switching. however, for type-i ifn responses a promoter activation in a cooperative manner has recently been discarded. taking the recent literature into account, we developed a model on the basis of a positive feedback loop circuit comprising the independent activation of isre and irfe elements by one isgf molecule and one irf dimer, respectively. this model has the needed and sufficient elements to sustain a bistable system. additionally, multiple model structures testing different biological scenarios including: additional feedback loops, additional intermediate elements, and promoter cooperative activation were tested, obtaining that the presented model reproduces the experimental data best. our simulation results agree with the experimental data obtained from flow cytometry. additional characteristics observed in the experimental data are also reproduced by the model, such as the basal expression of irf without ifn stimulation, see additional file the number of molecules in the active state of the system agrees with the order of magnitude reported for mammal mrnas and the range of protein concentration.  here, we present a method to fit stochastic models to experimental data. the method is based on the comparison of distributions. the central idea of the method is to use a deterministic precondition that is defined by the experimental data as a filter avoiding large amounts of costly stochastic simulations. using this idea, the number of parameters evaluated under stochastic dynamics is reduced, resulting in a significant improvement in the performance of the algorithm. as a case study, we used a model of irf gene expression investigating the origin of bimodality in its dynamics upon ifn stimulation. our results allowed us to conclude that a circuit with irf promoter activation by one irf dimer and one isgf molecule is sufficient to explain the observed bimodal dynamics. additional files additional file pseudocode for the genetic algorithm with the deterministic precondition. additional file with the pseudocode for the genetic algorithm. additional file identifying parameters for a constitutive gene expression circuit with in-silico data. additional file with an example where the method is applied to a simulation case with synthetic data. the file contains figures a1-a and tables a and a additional file random search strategy. additional file with figure a results obtained by using the random search strategy. additional file algorithm’s speed and performance analysis. additional file with an analysis of the algorithm’s speed and performance. the file contains figure a and tables a and a additional file testing the reproducibility of the selected parameter values. additional file containing figure a that shows the parameter distributions. additional file fitting the model to different doses of ifn. additional file with figure a and a that shows the fits of the model to flow cytometry data describing the stimulation of cells with 100u and 250u of ifn. additional file model dynamics without ifn stimulation. additional file with figure a that shows the time courses for the model in absence of ifn stimulation. abbreviations auarbitrary units bacbacterial artificial chromosome dpdeterministic precondition gagenetic algorithm ifninterferon isginterferon stimulated genes odeordinary differential equation pdfprobability density function rsrandom search sssstable steady state 
mathematical models are used to gain an integrative understanding of biochemical processes and networks. as a case study we fitted a model of the irf gene expression circuit to time-course experimental data obtained by flow cytometry. irf shows bimodal dynamics upon ifn stimulation. however, compared to the wealth of computational methods used to fit and analyze deterministic models, there is only little available to quantify the exactness of the fit of stochastic models compared to experimental data or to analyze different aspects of the modeling results. irf7ifnstochastic modelsparameter estimationconacyt263795bundesministerium für bildung und forschungimmunoquantissue-copyright-statement© the author 2017 commonly the models are based on deterministic ordinary differential equations. our results allow us to conclude that the activation of the irf promoter by the combination of irf and isgf is sufficient to explain the observed bimodal dynamics. when molecular counts are low, stochastic formalisms like monte carlo simulations are more appropriate and well established. the method uses a comparison of the probability density functions that are computed based on monte carlo simulations and the experimental data. the method improves its performance by selecting parameters values after comparing the similitude between the deterministic stability of the system and the modes in the experimental data distribution. here, we developed a method to fit stochastic simulations to experimental high-throughput data, meaning data that exhibits distributions.
 mathematical models are used to gain an integrative understanding of biochemical processes and networks. commonly the models are based on deterministic ordinary differential equations. when molecular counts are low, stochastic formalisms like monte carlo simulations are more appropriate and well established. however, compared to the wealth of computational methods used to fit and analyze deterministic models, there is only little available to quantify the exactness of the fit of stochastic models compared to experimental data or to analyze different aspects of the modeling results.  here, we developed a method to fit stochastic simulations to experimental high-throughput data, meaning data that exhibits distributions. the method uses a comparison of the probability density functions that are computed based on monte carlo simulations and the experimental data. multiple parameter values are iteratively evaluated using optimization routines. the method improves its performance by selecting parameters values after comparing the similitude between the deterministic stability of the system and the modes in the experimental data distribution. as a case study we fitted a model of the irf gene expression circuit to time-course experimental data obtained by flow cytometry. irf shows bimodal dynamics upon ifn stimulation. this dynamics occurs due to the switching between active and basal states of the irf promoter. however, the exact molecular mechanisms responsible for the bimodality of irf is not fully understood.  our results allow us to conclude that the activation of the irf promoter by the combination of irf and isgf is sufficient to explain the observed bimodal dynamics. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords irf7ifnstochastic modelsparameter estimationconacyt263795bundesministerium für bildung und forschungimmunoquantissue-copyright-statement© the author 2017 