Parameter Estimation Methods for Chaotic Intercellular Networks

We have investigated simulation-based techniques for parameter estimation in chaotic intercellular networks. The proposed methodology combines a synchronization–based framework for parameter estimation in coupled chaotic systems with some state–of–the–art computational inference methods borrowed from the field of computational statistics. The first method is a stochastic optimization algorithm, known as accelerated random search method, and the other two techniques are based on approximate Bayesian computation. The latter is a general methodology for non–parametric inference that can be applied to practically any system of interest. The first method based on approximate Bayesian computation is a Markov Chain Monte Carlo scheme that generates a series of random parameter realizations for which a low synchronization error is guaranteed. We show that accurate parameter estimates can be obtained by averaging over these realizations. The second ABC–based technique is a Sequential Monte Carlo scheme. The algorithm generates a sequence of “populations”, i.e., sets of randomly generated parameter values, where the members of a certain population attain a synchronization error that is lesser than the error attained by members of the previous population. Again, we show that accurate estimates can be obtained by averaging over the parameter values in the last population of the sequence. We have analysed how effective these methods are from a computational perspective. For the numerical simulations we have considered a network that consists of two modified repressilators with identical parameters, coupled by the fast diffusion of the autoinducer across the cell membranes.


Introduction
Most dynamical systems studied in the physical, biological and social sciences that exhibit a rich dynamical behavior can be modeled by sets of nonlinear differential equations. These mathematical models are a useful tool to predict complex behaviors using numerical simulations. However, for the vast majority of systems, and particularly for biological systems, we lack a reliable description of the parameters of the model. In this paper we are interested in parameter estimation for coupled intercellular networks displaying chaotic behavior, since models of spontaneously active neural circuits typically exhibit chaotic dynamics (for example, spiking models of spontaneous activity in cortical circuits [1][2][3] and the analogous spontaneously active firing-rate model networks [4,5]).
The problem of parameter estimation can be tackled in different ways, e.g., using multiple shooting methods [6][7][8] or some statistical procedures based on time discretizations and other approximations [9][10][11][12][13]. These methods involve the solution of high-dimensional minimization problems, since not only the unknown parameters but also the initial values of the trajectory segments between the sampling times need to be estimated [7,14]. This is specially difficult when working with chaotic systems since very complicated error landscapes with many local minima can appear. In particular, notice that chaotic systems have an exponential sensitivity to initial conditions, that is, completely different trajectories can be obtained for identical parameter values and very similar initial conditions of the system variables.
On the other hand, several authors have suggested to take advantage of synchronization techniques for coupled chaotic systems and turn them into accurate parameter estimation methods [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. There is a variety of techniques that rely on the synchronization properties of chaotic systems in order to tackle the parameter estimation problem. For example, some authors have proposed to handle the parameters as additional variables whose dynamics is described by tailored differential equations designed to have a fixed point at the true parameter values [14][15][16][17][18][19][20]. Adaptive estimation techniques relying on the time discretization of the state trajectories [22,23,30] and Monte Carlo methods [29,31], including particle filters [32,33], have also been investigated. All these techniques address the estimation of the system parameters independently of the initial conditions of the dynamic variables.
In this paper we investigate techniques that combine a synchronization-based framework for parameter estimation in coupled chaotic systems with some state-of-the-art computational inference methods borrowed from the recent literature in computational statistics. In particular, we describe estimation methods based on the accelerated random search (ARS) optimization algorithm [29,31,34] and two approximate Bayesian computation (ABC) [35][36][37][38][39] schemes. ABC is a general methodology for non-parametric inference that can be applied to practically any system of interest. We investigate two ABC-based methods. The first one is a Markov Chain Monte Carlo (MCMC) scheme [37] that generates a series of random parameter realizations for which a low synchronization error is guaranteed. We show that accurate parameter estimates can be obtained by averaging over these realizations. The second ABC scheme is termed Sequential Monte Carlo (SMC) ABC [38]. It generates a sequence of ''populations'', i.e., sets of randomly generated parameter values, where the members of the n{th population attain a synchronization error that is lesser than the error attained by members of the (n{1){th population. Again, we show how very accurate estimates can be obtained by averaging over the parameter values in the last population of the sequence. For the numerical simulations we consider a network of two coupled repressilators, since the repressilator is a paradigmatic gene regulatory system.

Structure of the Model
The repressilator is a prototype of a synthetic genetic clock built by three genes and the protein product of each gene represses the expression of another in a cyclic manner [40]. It can be constructed experimentally and produce near harmonic oscilla-tions in protein levels. In the original repressilator design [40], the gene lacI expresses protein LacI, which inhibits transcription of the gene tetR. The product of the latter, TetR, inhibits transcription of the gene cI. Finally, the protein product CI of the gene cI inhibits expression of lacI and completes the cycle. (See left-hand module in Fig. 1 of Ref. [41]).
Cell-to-cell communication was introduced to the repressilator by García-Ojalvo and coworkers [42] by introducing an additional feedback loop to the network scheme that is based on the quorum sensing mechanism. The additional genetic module, which might be placed on a separate plasmid involves two other proteins [42][43][44]: LuxI, which produces a small autoinducer (AI) molecule that can diffuse through the cell membrane, and LuxR, which responds to the autoinducer by activating transcription of a second copy of the repressilator gene lacI. The additional quorum sensing feedback loop can be connected to the basic repressilator in such a way that it reinforces the oscillations of the repressilator or competes with the overall negative feedback of the basic repressilator. The first one leads to phase attractive coupling for a robust synchronised oscillation [42] whereas the latter one evokes phase-repulsive influence [45][46][47], which is the key to multistability and a very rich dynamics including chaotic oscillations [41,48,49]. Placing the gene luxI under inhibitory control of the repressilator protein TetR (see Fig. 1 bottom in Ref. [41]) leads to the desired repressive and phase-repulsive coupling. We term this system modified repressilator. Phase repulsive coupling is common in several biological systems, e.g. in neural activity in the brain of songbirds [50], in the respiratory system [51], in the jamming avoidance response in electrical fish [52] and in the morphogenesis in Hydra regeneration and animal coat pattern formation [53].
In particular, in this work we are going to consider two modified repressilators with identical parameters, coupled by the fast diffusion of the autoinducer (AI) across the cell membranes. The mRNA dynamics is described by the following Hill-type kinetics with Hill coefficient n: where the subscript i~1,2 specifies the cell, and a i , b i , and c i represent the concentrations of mRNA molecules transcribed from the genes of tetR, cI, and lacI, respectively. The parameter a is the dimensionless transcription rate in the absence of a repressor. The parameter k is the maximum transcription rate of the LuxR promoter.
The protein dynamics is given by where variables A i , B i , and C i denote the concentration of the proteins TetR, CI, and LacI, respectively. The dynamics of the proteins is linked to the amount of the responsible mRNA, and the parameters b a,b,c describe the ratio between mRNA and the protein lifetimes (inverse degradation rates). In what follows, we are going to assume b b~bc . Thus, we can consider it as a single parameter. The model is made dimensionless by measuring time in units of the mRNA lifetime (assumed equal for all genes) and the mRNA and protein levels in units of their Michaelis constant. The mRNA concentrations are additionally rescaled by the ratio of their protein degradation and translation rates [42].
The third term on the right-hand side of Eq. (3) represents activated production of lacI by the AI molecule, whose concentration inside cell i is denoted by S i . The dynamics of CI and LuxI can be considered identical, assuming equal lifetimes of the two proteins and given that their production is controlled by the same protein (TetR). Hence, the synthesis of the AI S i can be considered to be controlled by the concentration B i of the protein CI. Taking also into account the intracellular degradation of the AI and its diffusion toward or from the intercellular space, the dynamics of S i is given by where the diffusion coefficient g depends on the permeability of the membrane to the AI. Because of the fast diffusion of the extracellular AI (S e ) compared to the repressilator period, we can apply the quasi-steady-state approximation to the dynamics of the external AI [42], which leads to The parameter Q is defined as where N~2 is the number of cells, V ext is the total extracellular volume, k se is the extracellular AI degradation rate, and d is the product of the membrane permeability and the surface area. We achieve chaotic behavior for the following parameter values [41]: Q~0:85, n~2:6, a~216, b a~0 :85, b b~0 :1, g~2, k~25, k s0~1 , and k s1~0 :01. In particular, these are the values we are considering throughout this manuscript in order to assess the parameter estimation algorithms.
To see how the system behaves around these values we have plotted some bifurcation diagrams taking in each one a different parameter as a control parameter, whereas the rest of the parameters remain constant on the values mentioned above. Figure 1(A-I) represents the bifurcation diagrams versus the different control parameters. In particular, in each plot we are representing all maxima of the variable a 1 for some fixed initial condition as a function of the corresponding control parameter. We can observe how for some values of the control parameters the behavior changes from periodic to chaotic and vice versa. The vertical red dashed lines correspond to the values we are considering to assess the parameter estimation algorithms. Notice how for these values the system can display chaotic behavior.

Problem Statement
We first introduce the notation to be used in the description of the parameter estimation methodologies. Let represent the chaotic intercellular network (consisting of two coupled modified repressilators), with state variables and unknown parameters p~(p 1 , . . . ,p 9 ) The vector-valued function f can be explicitly written down by comparing Eq. (10) with Eqs. (1)- (7). In the sequel we refer to the system of Eq. (10) as primary system. Since f in Eq. (10) is known, we can build a model with identical functional form but adjustable parameters and a coupling term, termed secondary system in the sequel, as where y~(y 1 , . . . ,y 14 ) is the time-varying vector that contains the model state variables, is the adjustable parameter vector, D is a coupling coefficient and s i,j : R 14 ?R 14 is a vector function that selects the i-th and the j-th element of its argument, i.e., The definition of the latter function implies that coupling is performed only through two scalar time series, x i and x j , from the primary system. In particular, the coupling scheme we have chosen for the simulation setup in this paper is _ a a a a i~{ a a i z a a 1z where i~1,2 denotes the cell number (i.e., the repressilator index). Thus, the coupling only appears in two of the fourteen differential equations we have in the secondary system. Since we assume identical functional form for the primary and secondary systems, when the secondary parameter vector, q, coincides with the primary parameter vector, p, the state variables y synchronize with x for DwD c , where D c is a coupling threshold [29]. On the contrary, if q=p then identical synchronization cannot occur. However, the difference y{x k k is expected to be small whenever the difference of the two parameter vectors is small and D is sufficiently large. Therefore, the objective of a parameter estimation method based on synchronization is to compute a valuê q q such that y{x k k&0 over time, since the latter impliesq q&p. Figure 2(left) shows how the synchronization error D a a 1 {a 1 D evolves in time when considering D~20 and identical parameters for the primary and the secondary systems. We can see how this error can be very small, less that 10 {10 . We are going to fix D~20 in all simulations throughout the paper. In Fig. 2(right) we are representing one variable of the secondary system (y 1~ a a 1 ) versus the same variable of the primary system (x 1~a1 ) and observing the typical straight line, characteristic of the synchronization phenomenon.
It has been verified by means of numerical simulations that in order to obtain synchronization at least two coupling variables are needed (excluding S 1 or S 2 ), one from each repressilator in the primary system. Therefore, other combinations different from (a 1 ,a 2 ) are also possible, as for example (b 1 ,A 2 ). Coupling via a single variable is not sufficient to guarantee that the secondary system synchronizes with the primary one.

Parameter Estimation Methods
Here we introduce three parameter estimation methods that combine the synchronization-based framework described above with state-of-the-art computational methods. In particular, we present results for the joint estimation of four parameters, in the coupled modified repressilator network. In all simulations shown in this manuscript we have assumed that just two scalar signals from the primary system are observed, namely the coupling variables (a 1 ,a 2 ), and we have numerically integrated the secondary system using a second order Runge-Kutta method with a time step T s~0 :005 time units (t.u.).

Accelerated Random Search
We first focus on a parameter estimation method for chaotic intercellular networks that takes advantage of chaos synchronization and is based on an efficient Monte Carlo optimization procedure, known as accelerated random search (ARS) method [29,31]. In particular, the method we are going to describe consists in a variation of the technique proposed in [31]. Assuming that the variables (a 1 ,a 2 ) are observed during the time interval (0,T o ), the value of the parameters in the secondary system can be calculated as the solution to the optimization problem where q~(q 1 , . . . ,q 4 )~( Q Q, b b a , n n, a a) and the cost function is a quantitative representation of the synchronization error between the primary and the secondary systems. Notice that N~T o =T s , where T o~1 0 5 t.u. The method can be outlined as follows: 1. Initialization. Choose: (a) initial parameter valueŝ q q(0)~(q q 1 (0),q q 2 (0),:: where m represents the number of parameters we want to estimate; (b) maximum and minumum ''search variances'' for each parameter, s 2 j,max and s 2 j,min , respectively, where j~1, . . . ,m; (c) a ''contraction factor'' for each parameter, c j w1, j~1, . . . ,m.
The algorithm can be stopped when the synchronization error J(q q(n)) reaches a certain threshold or after a prescribed number of steps n (e.g., when n~n o for some sufficiently large n o ). Notice that the total number of iterations is n o L and the time evolution of the secondary system state, y, has to be integrated at each iteration, since a new candidate parameter vector is drawn each time.
We have carried out numerical simulations where this Monte Carlo optimization algorithm has been iterated a total of 8|10 4 times, with the number of consecutive iterations for each parameter (iterations per loop) set to L~100, and the secondary system has been integrated for each iteration during T o~1 The normalized quadratic error for each parameter q i , i~1, . . . ,4, is defined as,

Approximate Bayesian Computation
In Bayesian estimation theory the parameters are modeled as random variables, rather than unknown but deterministic numbers. Consequently, ABC-based methods aim at approximating the probability distribution of the parameters, vector q, conditional on the observations from the primary system. The technique, to be described next, is nonparametric, i.e., the distribution of q is approximated by a set of random samples.
ABC methods have been conceived with the aim of inferring posterior distributions without having to compute likelihood functions [35][36][37][38][39]. The calculation of the likelihood is replaced by a comparison between the observed and the simulated data. In the setup of this paper, the comparison is carried out between the data from the primary system (the observations) and the data generated by the secondary system (the model with adjustable parameters). The comparison between these data represents, in our case, a measure of the synchronization error between these two systems.
Let p denote the a priori probability density function (pdf) of the random parameter vector q, let f (yDq) denote the probability distribution of the data y generated by the secondary system conditional on the realization of q and let d(x,y) be a distance function that measures the synchronization error by comparing the observed time series x from the primary system and the generated time series y from the secondary one. Since the system of interest is deterministic, f (yDq) collapses into a Dirac delta measure when q is given. In the ABC framework, though, we are not interested in evaluating f but rather in generating y using the model conditional on q. This amounts to integrating the secondary system with a given realization of the adjustable parameters q.
The simplest ABC algorithm is the ABC rejection sampler [35]. For a given synchronization error threshold (often termed tolerance) w0, the algorithm can be described as follows. 1. Sampleq q from p(q). 2. Simulate a datasetỹ y from the conditional probability distribution f (yDq q).

If the distance function (synchronization error) is d(x,ỹ y)ƒE,
then acceptq q, otherwise reject it. 4. Return to the first step.
The output of an ABC algorithm is a population of parameter values randomly drawn from the distribution p(qDd(x,ỹ y)ƒ), i.e., a distribution with density proportional to the prior but restricted to the set of values of q for which the synchronization error is at most .
The disadvantage of the ABC rejection sampler is that the acceptance rate is very low when the set of values of q for which the synchronization error is less than turns out to be relatively small. Thus, instead of implementing this method we have decided to implement two more sophisticated ABC algorithms. The first one is a Markov chain Monte Carlo (ABC MCMC) technique and the second one is a sequential Monte Carlo (ABC SMC) method.

Approximate Bayesian Computation Markov Chain Monte Carlo
The Approximate Bayesian Computation Markov Chain Monte Carlo (ABC MCMC) algorithm is a Metropolis-Hastings [37] MCMC method that incorporates one additional test to ensure that all parameters in the chain yield a synchronization error below the threshold . It can be outlined as follows.
The outcome of this algorithm is a Markov chain with the stationary distribution p(qDd(x,ỹ y)ƒE) [37]. The parameters are assumed independent a priori, hence p(q)~P Since the proposal distribution is symmetric, g(q i Dq q)~g(q qDq i ), and the prior is uniform, the acceptance probability in Eq.(32) is j~1.
The distance function (synchronization error) has been chosen as where T o~1 0 5 t.u. and NT s~T0 . This expression is equivalent to the cost function in the ARS method. The tolerance (threshold for the synchronization error) has been chosen as E~0:001. A chain of 5|10 5 samples has been generated, what implies that the secondary system has been integrated 5|10 5 times. The initial point of the chain is selected to ensure that the associated distance is less than 0:05. Figure 4(A-D) shows the histograms of the approximate marginal posterior distributions for each parameter. In order to reduce the strong correlation between consecutive samples in the Markov chain we have subsampled by a factor of 50. We have calculated the mean values of each histogram as well as the normalized quadratic errors according to the following expression where q q i , for i~1, . . . ,4, represents the mean value of the histogram of the corresponding parameter. The values of the normalized quadratic errors for the estimated parameters are We can see how three parameters are accurately estimated whereas for one of them, b a the error is significantly higher compared to (31).

Approximate Bayesian Computation Sequential Monte Carlo
A more sophisticated application of the ABC methodology is the Approximate Bayesian Computation Sequential Monte Carlo algorithm [39,54,55]. In ABC SMC, a number of sampled parameter values (often termed particles), fq (1) , . . . ,q (N) g, drawn from the prior distribution p(q), are propagated through a sequence of intermediate distributions p(qDd(x,ỹ y)ƒ i ), i~1, . . . ,T{1, until they are converted into samples from the target distribution p(qDd(x,ỹ y)ƒ T ). The tolerances are chosen such that 1 w . . . w T w0, thus the empirical distributions gradually evolve towards the target posterior. The ABC SMC algorithm proceeds as follows [39].
The prior distributions we have considered for each parameter are p(Q)~p(b a )~p(n)~U(0,5) and p(a)~U(0,500), the same as for the ABC MCMC algorithm. The perturbation kernel K t is Gaussian, namely  Figure 5(A-D) shows the histograms of the approximate marginal posterior distributions for each parameter where the different weights of the different particles have been taken into account for the representation. We have calculated the mean values of each histogram (that match the actual values) as well as the normalized quadratic errors using Eq. (38), that is, the same expression as for the ABC MCMC algorithm. The values of the normalized quadratic errors for the estimated parameters with the ABC SMC method are Figure 6(A-D) shows the output (i.e. the accepted particles) of the ABC SMC algorithm as scatterplots of some of the twodimensional parameter combinations, where we have information of different populations in the same plot. As we iterate the algorithm we obtain populations that are more dense around the desired values, as shown in these plots, where the particles from the prior are represented in blue, particles from population 2 in green, particles from population 4 in light blue, population 6 in  In order to gain insight of how the parameters are estimated during the evolution of the algorithm, we have represented some box-plot diagrams, one for each parameter to be estimated, as seen in fig. 7(A-D). In each diagram, we have information about the corresponding parameter as a function of the population index. In particular, the central mark of each box is the median of the population, the edges of the box are the 25-th and 75-th percentiles, the whiskers extend to the most extreme data points not considered as outliers, and the outliers are plotted individually using the plus symbols in red. The horizontal lines in red represent the actual values of the parameters we have used in our simulation. We can see from these plots how for a high enough population index the median values of the four parameters perfectly match the actual values.
We have also studied the computational cost of the algorithm. In fig. 8 we can see the number of samples or particles we have generated in order to have N~400 accepted particles for each population. We can see how the number of particles increases with the population index, being significantly high for the last population index, since it corresponds to a very small value of the synchronization error. Notice that the vertical axis of this figure is in a logarithmic scale.

Comparison of the Methods
Here we compare not only the accuracy but also the computational complexity of all three Monte Carlo methods for the joint estimation of the four parameters, Q, b a , n and a, of the chaotic intercellular network. To do that, we have calculated for the ABC SMC algorithm the computational load up to each population. Specifically, the computational complexity of generating a sequence of m populations is given by the number of samples or particles that have to be generated before completing the m{th population. Note that the computational load for the m{th population also includes all samples needed to generate the m{1 previous populations. Fig. 9 provides a graphical depiction of the complexity of the three methods, that we have investigated (ARS, ABC MCMC and ABC SMC algorithms). The line in red represents the complexity for the ARS method and the line in blue indicates the complexity for the ABC MCMC technique.
The parameter estimation errors attained with the ARS method are of the same order as the errors of the parameter estimates computed from the 12{th population of the ABC SMC algorithm. However, the number of samples generated to run the ARS procedure is &3|10 4 while the ABC SMC technique demands the generation of &4|10 5 random samples up to the 12{th population. The ABC MCMC method achieves the poorest performance, as it requires the generation of the highest number of samples (&5|10 5 ) and produces the largest errors (up to three orders of magnitude worse than the ARS or ABC SMC estimates).

Discussion
We have investigated three computational inference techniques for parameter estimation in a chaotic intercellular network that consists of two coupled modified repressilators. The proposed methodology combines a synchronization-based framework for parameter estimation in coupled chaotic systems with some stateof-the-art computational inference methods borrowed from computational statistics. In particular, we have focussed on an accelerated random search algorithm and two approximate Bayesian computation schemes (ABC MCMC and ABC SMC). The three methods exploit the synchronization property of chaotic systems. Therefore, it is not necessary to estimate the initial conditions of the variables, which is an important advantage from a computational point of view.
We have carried out the numerical study in this paper assuming that only two variables from the primary system can be observed. This is the minimum number of observed variables in order to guarantee the synchronization of the secondary system. If additional variables can be observed it is possible to easily incorporate them into the proposal methodology. For example, if the variables x 1 , . . . ,x 8 are observed (this is the full state of the first repressilator and the first variable, x 8 , of the second repressilator) we can redefine the distance function of Eq. (38) as It can be verified (numerically) that using the distance in (44) (which intuitively provides ''more information'' about the primary system) leads to more accurate parameter estimates (or, alternatively, a greater number of parameters can be estimated if necessary). Note, however, that this comes at the expense of an additional computational effort and, moreover, it is unclear that all these variables can be accurately measured in practice.
The proposed methods can be applied when the observed time series are contaminated with additive noise of moderate variance. For example, if the observations have the formã a 1 (nT s )~a 1 (nT s )zu 1 (n) andã a 2 (nT s )~a 2 (nT s )zu 2 (n), where u 1 (n) and u 2 (n) are sequences of independent and identically distributed Gaussian noise variables with zero mean and variance s 2 u , then the distance function of Eq. (30) is lower-bounded by the noise variance. Specifically, ifq q&p and, hence, we assume that a a i (nT s )&a i (nT s ), it turns out that the distance d(y,x) is an estimator of (twice) the noise variance ( a a 1 (nT s ){ã a 1 (nT s )) 2 z( a a 2 (nT s ){ã a 2 (nT s )) 2 À Á This indicates that the synchronization error cannot go below the (approximate) bound of 2s 2 u and, therefore, the ABC-based methods can work as long as the tolerances ( t , t~1, . . . ,T, in the ABC SMC method, or E in the ABC MCMC technique) are chosen to be greater than 2s 2 u . This means that the ABC SMC algorithm with T~12 populations and T~1 0 {4 can still provide accurate parameter estimates when the observation noise variance is s 2 u v T =2. In order to handle larger noise variances, one needs to relax the coupling (i.e., choose a smaller coupling factor D) and increase the observation period T o . This makes the distance function more sensitive to the discrepancy between q and p, which in practice means that we can choose a larger tolerance (e.g, T in the ABC SMC algorithm) and preserve the accuracy of the resulting estimateq q.