Automated parameter estimation for biological models using Bayesian statistical model checking

Parameter estimation using global and local search

Sun et al. [18] survey metaheuristic techniques used in parameter estimation in systems biology, focusing on simulated annealing, evolutionary algorithms and hybrid strategies that combine multiple heuristics.

Gonzalez et al. [19] use simulated annealing [20] to find parameters in S-system models of biochemical networks. At a given parameter point during the annealing process, they find a neighboring point by adding a noise term to each component of the parameter vector that is dependent on the current optimization error.

Lillacci and Khammash have designed a method that uses Kalman filtering, statistical testing and numerical optimization for parameter estimation and model selection [21]. Inspired by these two approaches [19, 21], we have used simulated annealing for optimization and statistical hypothesis testing-based verification in our parameter estimation algorithm for probabilistic models.

Algorithms based on maximum-likelihood estimation (MLE) and the singular value decomposition (SVD) have been proposed by Reinker et al. [22], to estimate the reaction rate constants (the parameters) from discrete time series data for molecule counts in stochastic biochemical reactions.

Rodriquez-Fernandez et al. use a hybrid approach for parameter estimation in deterministic, non-linear models of biochemical pathways [23]. Their approach combining local and global optimization methods helps overcome the problem of convergence to local minima common in traditional local optimization methods (like gradient-descent) and the problem of slow convergence seen in global optimization techniques.

Moles et al. have studied the global optimization based approaches to the parameter estimation problem in nonlinear dynamic systems using a 36 parameter dynamic pathway model as a benchmark [24]. They report that only approaches based on evolutionary strategies, viz. the stochastic ranking evolutionary strategy (SRES) and unconstrained evolutionary strategy (use), were able to successfully find parameters in the pathway.

The parameter estimation problem for biological pathway models has also been addressed by Koh et al. [25]. An interesting element of their approach is to decompose the model into components whose parameters can be estimated independently. They model pathways using hybrid Petri nets and use evolutionary strategies for parameter estimation, but their approach can be used for any modeling framework and choice of the optimization technique.

Matthew et al. have used global sensitivity analysis (GSA) to study the effect of parameter perturbations in a large, non-linear, ODE-based model of acute inflammation in mice [26], and found that Interleukin 6 (IL-6) and nitric oxide (NO) had a significant, non-linear impact on inflammatory damage [27].

In another approach that uses sensitivity analysis, Donze et al. have developed an algorithm for parameter synthesis in nonlinear systems of ODEs (ordinary differential equations) [28], and applied their technique to find parameters in two models of acute inflammation [29, 30].

Next, we discuss parameter estimation approaches based on formal verification techniques.

Parameter estimation using formal verification

Many recent procedures for parameter estimation make use of a formal verification technique known as model checking.

A model checking approach to finding parameters in biological models has been used by Calzone et al [31]. They use the Biochemical Abstract Machine (BIOCHAM) modeling framework to describe the system and temporal logic model checking [32] to find parameters values in a user-specified range.

Dreossi and Dang have recently designed an algorithm that reduces the problem of parameter synthesis in polynomial (discrete-time) dynamical systems [33] to solving a set of linear programs [34]. They use their technique to find parameters of two well-studied epidemic models.

Batt et al. have used symbolic model checking to find parameters in a piecewise affine differential equation (PADE) model of the gene regulation IRMA network [35]. IRMA stands for in vivo "benchmarking" of reverse-engineering and modeling approaches [36]. Donaldson and Gilbert have designed a technique for parameter estimation that combines model checking with a genetic algorithm [37], and used it for estimating kinetic rate constants in a model of the mitogenactivated protein kinases (MAPK) signaling pathway.

Usually, model checking based methods of parameter estimation require that the relevant specification be expressed in a formal temporal logic, whose satisfaction against a given execution path (known as a trace) of the model can be determined. Rizk et al. define a continuous degree of satisfaction of a temporal logic formula for any given trace of the model [38] and use it as a fitness function to drive an optimization-style search for kinetic parameters in models of the cell cycle and the MAPK signal transduction.

In recent work, Mancini et al. have used statistical model checking to synthesize parameters in an ODE-based biological model [39]. They have implemented a distributed, multi-core version of their algorithmic technique, and used it to estimate patient-specific parameters of a a human menstrual cycle model [40].

Stochastic biological models

Many biological systems have traditionally been described using deterministic, mathematical models, often in terms of ordinary differential equations [[12] §2.1]. However, in recent years, the need for incorporating the uncertainty inherent in biological systems has led to the development of stochastic models: these are harder to analyze but more accurately reflect the behavior of the underlying system [5]. Common stochastic models in biology include Markov jump processes [5], stochastic differential equations [41], discrete-time Markov chains [[42], Chapter 3] and continuous-time Markov chains [43].

Also, researchers are increasingly developing computational models that naturally capture the bottom-up nature of biological phenomena and hence are more amenable to in-silico implementation [44]. Such models are constructed by observing large sets of timeseries data, combined with their expert insight into the system being modeled [45]. Some of these models are based on experimental data of varying veracity and error propagation into the designed model is an ever-present challenge [46].

We will focus on discrete-time Markov chains, a class of stochastic models that are widely used in the sciences, engineering, economics and other areas. We closely follow Baier and Katoen [8] in formally defining them.

Definition 1 (Discrete-time Markov chain) A discrete time Markov chain (DTMC) is given by $\mathcal{M}=\left(S,P,init,AP,L\right)$ where:

We next define (a) a parameterized family of DTMCs, given a set of parameters and (b) execution paths over them.

Definition 2 (Parameterized DTMC) A parameterized discrete time Markov chain (ParDTMC) is given by $\mathcal{M}=\left(S,{\Theta }_m,P,init,AP,L\right)$ where

Definition 3 (ParDTMC execution paths) An execution path of a ParDTMC $\mathcal{M}=\left(S,{\Theta }_m,P,init,AP,L\right)$ given by σ = s₀, s₁, . . .. The suffix of a path σ that starts at state s _i (i.e. s _i , s_i+1, . . .) is denoted σ ⁱ .

As a case study, we have applied our automated parameter estimation technique to an agent-based model of the acute inflammatory response. We now briefly discuss major characteristics of ABMs.

Agent-based modeling in systems biology

Agent-based models (ABMs) form an interesting, wellstudied subset of complex probabilistic models. In recent years, agent-based modeling has emerged as a popular method for the representation, analysis and simulation of biological models [10].

ABMs allow the specification of high-level model properties as well as fine-grained component behavior. An ABM is composed of autonomous elements whose individual properties can be specified. At a macro-level, model-wide agent-interaction rules can be defined and enforced. Each agent has a physical location and its state evolves with time based on messages exchanged with other agents, allowing rich spatiotemporal properties to emerge [10].

ABMs can incorporate parallelism, object-oriented behavior and stochasticity, making them conducive to the development of computational models in systems biology. Another advantage of using ABMs is that they are bottom-up models that mirror the natural behavior of biological systems [47]. Agents interact with one another based on fixed rules and also have a spatial location [10].

Most agent-based models do not have compact analytic descriptions, and simulations across the input space must be performed to in order to infer general model properties [48]. Also, agent-based models tend to have a large number of parameters, some of them highly-sensitive in the sense that a small change can radically alter model behavior [49]. To the best of our knowledge, the parameter estimation problem in ABMs has only been addressed for very simple models [50] and most approaches use standard optimization techniques [31, 51].

As a case study for our parameter estimation technique, we have used an agent-based model of the acute inflammatory response due to endotoxin administration [52] written using the SPARK tool [53]. SPARK is an ABM framework designed for multi-scale modeling of biomedical models that is implemented in Java and allows users to develop models using its own programming language (SPARK-PL) [54–57].

One we have a basic model design available, we need to find the exact model instantiation in terms of concrete parameter values that meets desired behavior (usually experimental data). Therefore, we need a way to check if the models all requirements. We now discuss how to formally specify and automatically verify probabilistic computational models.

An ABM $\mathcal{A}$ consists of a fixed number of agents A₁, . . . A _n , where the state of agent A _i is determined by the values of variables $V_1^i\cdots V_m^i$. Assuming that for any agent A _i , i ∈ {1 . . . n}, variable $V_j^i$, j ∈ {1 . . . m} can take values in the set V ar _j , $\mathcal{A}$ can be represented as a DTMC.

Definition 4 (DTMC corresponding to an ABM) Consider ABM $\mathcal{A}\left(n,m\right)$ with n agents A₁, . . . A _n and m variables V₁ . . . V _m that take values over countably finite sets V ar₁, . . ., V ar _m respectively. The DTMC ${\mathcal{M}}_{\mathcal{A}}=\left(S,P,init,AP,L\right)$ corresponding to $\mathcal{A}\left(n,m\right)$ is:

We generalize agent-variable ABMs, by adding the notion of parameters to obtain a family of ABMs. More formally, we talk of a parametric ABM $\mathcal{A}\left(n,m,k\right)$ over n agents, m variables and k parameters whose dynamic behavior (i.e. its transition function P) now also depends on the parameter values. For a given parametric ABM, we can defined an equivalent ParDTMC.

Definition 5 (ParDTMC corresponding to a parametric ABM) Consider ABM $\mathcal{A}\left(n,m,k\right)$ with n agents A₁, . . . A _n and m variables V₁ . . . V _m that take values over countably finite sets V ar₁ . . . V ar _m and k parameters θ₁ . . . θ _k . The ParDTMC ${\mathcal{M}}_{\mathcal{A}}=\left(S,{\Theta }_k,P,init,AP,L\right)$ corresponding to $\mathcal{A}\left(n,m,k\right)$ is:

Remark 2 In order to use statistical model checking to solve the probabilistic model checking problem, we need to associate probabilities with executions paths of the model. For DTMCs, this has been shown by Kwiatkowska et al [9] and Baier and Katoen [[8], Chapter 10]. Although we do not prove this, this is also true for ParDTMCs. From now on, we assume that our models can be represented as ParDTMCs and thus have a unique underlying probability measure.

Specifying and checking biological properties

Our algorithm discovers parameters of stochastic biological from experimental data and expert behavioral specifications. While exploring the parameter space, we continually verify whether we have a parameter assignment at which the model meets the given specifications.

Typically, properties of biological systems that need to be formally specified are qualitative, quantitative and time-varying. Since the models are usually stochastic, specifications too should allow reasonable deviations from desired behavior. To capture such rich behavioral properties, we use a probabilistic temporal logic (see Definition 8) to specify expected model behavior.

In order to develop automated techniques for analyzing biological models, we also need the properties to be specified in a way that can be monitored as we perform ABM simulations. Monitoring is the process of determining if an execution trace of the model satisfies given specifications [58, 59]. Therefore, we write specifications using a logic belonging to a class of languages for which monitors can be derived algorithmically [60].

We translate natural language expert insights representing desired model behavior into a logic whose sentences are adapted finitely monitorable (AFM) formulas. The truth value of an AFM formula with respect to a trace of a model execution can be determined by observing a finite prefix of the trace [[61], §2],[[60], §2.2]. Next, we formally define our specification language:

Definition 6 (BLTL grammar) Given a ParDTMC $\mathcal{M}=\left(S,{\Theta }_k,P,init,AP,L\right)$, the grammar of bounded linear temporal logic (BLTL) formulas φ in Backus-Naur Form is as follows:

$⟨\varphi ⟩::=a|\varphi \wedge \varphi |\varphi \vee \varphi |\neg \varphi |\varphi {\text{U}}^d\varphi$

where a ∈ AP, $d\in \mathbb{N}$, and ∧, ∨ and ¬ are the usual proposition logic operators. We call U ^d the bounded until operator. ■

The semantics of a BLTL formula φ is defined over an execution path of a ParDTMC $\mathcal{M}$ (see Definition 3).

Definition 7 (BLTL semantics) Given a path σ of a ParDTMC $\mathcal{M}=\left(S,P,init,AP,L\right)$, the satisfaction of a BLTL formula φ on a path suffix ${\sigma }^i=\left(s_1,s_2,\dots \right)$ is denoted ${\sigma }^i\models \varphi$, and is determined by the following rules:

We now state, without proof, the fact that it is possible to check if any path of a ParDTMC satisfies a given BLTL formula, by observing a finite prefix of the path. For a proof of a similar lemma for continuous-time Markov chains, see Legay et al [62].

Lemma 1 (Monitorability of BLTL formulas over DTMCs) It is always possible to algorithmically decide if any simulation path σ = (s₁, s₂, . . .) of a ParDTMC $\mathcal{M}$satisfies given BLTL formula φ by observing a finite prefix of σ.

Note 1 (Additional bounded operators) We can define bounded versions of the usual linear temporal logic operators G (always) and F (eventually) as follows: F ^d ψ = true U ^d ψ, G ^d ψ = ¬F ^d ¬ψ. F ^d ψ means that ψ holds at some state within the next d state-transitions, and G ^d ψ means that ψ continually holds for the next d state-transitions.

The expressive power of our agent-based models emanates from their ability to capture uncertainty. For such stochastic models, a single execution trace that violates a given property cannot serve as a counterexample. We need a more flexible specification language whose formulas allow reasonable deviations from expected model behavior. For this, we define a probabilistic version of BLTL.

Definition 8 (Probabilistic Bounded Linear Temporal Logic) A specification of the form P _{≥ θ} (φ) is a probabilistic bounded linear temporal logic (PBLTL) formula if φ is a bounded linear temporal logic formula and θ is a probability threshold such that $\theta \in ℝ$, 0 ≤ θ ≤ 1. If a θ fraction of the traces of a model satisfy φ, we deem the model to satisfy the (probabilistic) AFM specification P _{≥ θ} (φ).■

Automated verification using model checking

Model checking is an automated technique for verifying finite and infinite state transition systems that is widely used for formal assurance of safety-critical systems [8]. Techniques based on model checking have been used for verification in a number of areas such as hybrid dynamical systems [63], computer software [64] and systems biology [43].

In order to use model checking to verify a time-varying system, the model is described using a Kripke structure and the property to be checked is written in a formal temporal logic [32]. We will use PBLTL (Definition 8) to define the probabilistic model checking problem.

Definition 9 (Probabilistic model checking (PMC)) Given a probabilistic model $\mathcal{M}$ and a PBLTL specification P _{≥ θ} (φ), determine if $\mathcal{M}$ satisfies φ with probability at least θ. ■

There are two approaches for solving the PMC problem:

Statistical approaches for probabilistic model checking are more scalable since they avoid the expensive calculation associated with accurately estimating exact probabilities. However, a limitation of statistical model checking algorithms is that the reported result is not guaranteed to be correct, i.e., there may be false positives or false negatives [71]. For an excellent overview of major statistical model checking techniques, we refer the reader to Legay et al [61].

Since computational models in systems biology often have large parameter spaces, we use statistical model checking as part of our parameter estimation algorithm. Next, we describe a reformulation of the probabilistic model checking problem in terms of statistical hypothesis testing - an approach first used by Younes [71].

Hypothesis testing based statistical model checking

We are interested in determining an answer to the probabilistic model checking problem, i.e. "Does $\mathcal{M}\models P_{\ge \theta }\left(\varphi \right)$?" (see Definition 9).

Assuming that the model's actual probability of satisfying the specification is u ∈[0, 1], i.e. $\mathcal{M}\models P_{=u}\left(\varphi \right)$, we test the (null) hypothesis H : u ≥ θ against the (alternative) hypothesis K : u <θ. If K is rejected we conclude that $\mathcal{M}\models P_{\ge \theta }\left(\varphi \right)$. If H is rejected we conclude that $\mathcal{M}⊭P_{\ge \theta }\left(\varphi \right)$.

Note 2 (Errors in hypothesis testing) For a hypothesis test procedure, the critical region is the part of the sample space where the null hypothesis is rejected. If the test rejects H when its true, it is considered a type-I error. On the other hand, if the test accepts H when its false, it is a type-2 error. Once the critical region for a test procedure has been decided, it uniquely determines the probabilities of type-1 and type-2 errors. For a given critical region, we denote by α (resp. β) the probability of a type-1 (resp. type-2) error [[72], §1.3].■

Naturally, for any statistical hypothesis testing procedure we want a critical region that minimizes the probabilities α and β of Type-1 and Type-2 errors (see Note 2). However, this implies either using a large value for either α or β, or drawing a large number of samples to ensure test accuracy.

Younes [70] suggests that the solution is to use the more relaxed test of H₀ : u ≥ u _r against H₁ : u ≤ u _l , where 0 ≤ u _l <θ <u _r ≤ 1. If H₀ (H₁) is accepted, we consider H (resp. K) to be true.

Remark 3 (Indifierence regions) [u _l , u _r ] is known as the indifierence region. If u ∈ [u _l , u _r ] (i.e. when both H₀ and H₁ are false), we do not care about the result of the test; thus the test procedure is allowed to accept either hypothesis. In practice, we often choose the indifierence region [u _l , u _r ] to be of width 2 * ∈ (where $0<\in \ll 1$), i.e. [u - ∈, u + ∈].

Our parameter estimation technique uses the Bayesian statistical model checking (BSMC) algorithm developed by Jha at al [16], to algorithmically check if a probabilistic model satisfies given behavioral specifications. Note that a fully Bayesian statistical model checking algorithm that verifies Dynamic Bayesian Networks against probabilistic queries was developed by Langmead [73].

Below, we briefly describe the main idea behind BSMC and refer to the reader to the literature for details [60, 61, 16].

Bayesian statistical model checking

Recall that we are interested in testing whether a probabilistic model meets a PBLTL specification with a minimum threshold probability, i.e. "Does $\mathcal{M}\models P_{\ge \theta }\left(\varphi \right)$?"

Assuming that $\mathcal{M}\models P_{=u}\left(\varphi \right)$, we have posed this problem as hypothesis testing query: test H : u ≥ θ against K : u <θ, and then relaxed it to use indifference regions: H₀ : u ≥ u _r against H₁ : u ≤ u _l (where 0 ≤ u _l < θ < u _r ≤ 1).

If H₁ is rejected, we consider that H : u ≥ θ holds, and hence conclude that $\mathcal{M}\models P_{\ge \theta }\left(\varphi \right)$. If H₀ is rejected, we conclude that $\mathcal{M}⊭P_{\ge \theta }\left(\varphi \right)$, i.e. $\mathcal{M}\models P_{<\theta }\left(\varphi \right)$.

Recall that we do not know the value of the actual probability u with which the model satisfies the specification i.e. $\mathcal{M}\models P_{=u}\left(\varphi \right)$. For Bayesian testing, we model this unknown probability as a random variable U and assume the availability of a prior density function g (.) that incorporates our existing knowledge of U.

The Bayesian statistical model checking procedure (see Algorithm 1) sequentially draws traces from $\mathcal{M}$ in an i.i.d. fashion until it rejects either H₀ or H₁. For each sample trace σ _i , (i ∈ {0, 1, . . .}) of $\mathcal{M}$, it uses a monitoring algorithm to determine if the path satisfies the BLTL specification φ.

Algorithm 1 Bayesian statistical model checking (BSMC) algorithm for probabilistic models

Require:

Probabilistic model $\mathcal{M}$,

PBLTL specification $P_{\ge \theta }\left(\varphi \right)$,

Threshold L > 1,

Prior density function g(·) of the unknown parameter u where $\mathcal{M}\models P_{=u}\left(\varphi \right)$

Indifference region bounds: (∈₁, ∈₂), where ∈₁ > 0, ∈₂ > 0.

{Note: Indifference region is [θ - ∈₁, θ + ∈₂]

{Note: H₀ : u ≥ θ + ∈₂; H₁ : u ≤ θ - ∈₁}

Ensure:

ans = false if H₀ rejected,

ans = true if H₀ is accepted,

n = Total number of traces sampled from $\mathcal{M}$.

1: n ← 0 {Number of traces drawn from $\mathcal{M}$.}

2: z ← 0 {Number of traces satisfying φ.}

3: repeat

4: Draw sample σ from $\mathcal{M}$

5: n ← n + 1

6: if σ ⊨ φ then

7: z ← z + 1

8: end if

9: $B\leftarrow \frac{{\int }_{\theta +{\in }_2}^1{u}^z{\left(1-u\right)}^{n-z}g\left(u\right)du}{{\int }_0^{\theta -{\in }_1}{u}^z{\left(1-u\right)}^{n-z}g\left(u\right)du}$ (2)

10: until $\left(B>L\right)\vee \left(B<\frac{1}{L}\right)$

11: if B >L then

12: ans ← true

13: else

14: ans ← f alse

15: end if

16: return (ans; n)

In summary, calculating the Bayes factor requires knowing the number of total traces n drawn from $\mathcal{M}$, the number z of traces that satisfied specification φ, the indifference region [u _l , u _r ] and the prior density g(·) of the unknown probability u. For more details about the Bayesian model checking algorithm, we refer the reader to the work of Jha [[60], Chapter 4].