Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 7-Day Trial for You or Your Team.

Learn More →

Exploring behaviors of stochastic differential equation models of biological systems using change of measures

Exploring behaviors of stochastic differential equation models of biological systems using change... Stochastic Differential Equations (SDE) are often used to model the stochastic dynamics of biological systems. Unfortunately, rare but biologically interesting behaviors (e.g., oncogenesis) can be difficult to observe in stochastic models. Consequently, the analysis of behaviors of SDE models using numerical simulations can be challenging. We introduce a method for solving the following problem: given a SDE model and a high-level behavioral specification about the dynamics of the model, algorithmically decide whether the model satisfies the specification. While there are a number of techniques for addressing this problem for discrete-state stochastic models, the analysis of SDE and other continuous-state models has received less attention. Our proposed solution uses a combination of Bayesian sequential hypothesis testing, non-identically distributed samples, and Girsanov’s theorem for change of measures to examine rare behaviors. We use our algorithm to analyze two SDE models of tumor dynamics. Our use of non-identically distributed samples sampling contributes to the state of the art in statistical verification and model checking of stochastic models by providing an effective means for exposing rare events in SDEs, while retaining the ability to compute bounds on the probability that those events occur. Background level description of a dynamical behavior, j (e.g., a for- The dynamics of biological systems are largely driven by mula in temporal logic). It then computes a statistically stochastic processes and subject to random external per- rigorous bound on the probability that the given model turbations. The consequences of such random processes exhibits the stated behavior using a combination of are often investigated through the development and ana- biased sampling and Bayesian Statistical Model Check- lysis of stochastic models (e.g., [1-4]). Unfortunately, the ing [7,8]. validation and analysis of stochastic models can be very Existing methods for validating and analyzing stochas- challenging [5,6], especially when the model is intended tic models often require extensive Monte Carlo sam- to investigate rare, but biologically significant behaviors pling of independent trajectories to verify that the (e.g., oncogenesis). The goal of this paper is to introduce model is consistent with known data, and to character- an algorithm for examining such rare behaviors in Sto- ize the model’s expected behavior under various initial chastic Differential Equation (SDE) models. The algo- conditions. Sampling strategies are either unbiased or rithm takes as input the SDE model, ,and ahigh- biased. Unbiased sampling strategies draw trajectories according to the probability distribution implied by the model, and are thus not well-suited to investigating rare * Correspondence: [email protected]; [email protected] 1 behaviors. For example, if the actual probability that the Electrical Engineering and Computer Science Department, University of -10 Central Florida, Orlando FL 32816 USA model will exhibit a given behavior is 10 , then the Computer Science Department, Carnegie Mellon University, Pittsburgh PA expected number of samples need to see such behaviors 15213 USA is about 10 (See Figure 1). Biased sampling strategies, Full list of author information is available at the end of the article © 2012 Jha and Langmead; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 2 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 Figure 1 Observing rare behaviors in i.i.d. sampling is challenging. A toy model with a one low-probability state. An unbiased sampling algorithm may require billions of samples in order to observe the ‘bad’ state. Statistical algorithms based on i.i.d. algorithms are not suitable for analyzing such models with rare interesting behaviors. in contrast, can be used to increase the probability of [7,8,20-22]). There are two key differences between the observing rare events, at the expense of distorting the method in this paper and existing methods. First, exist- underlying probability distribution/measure. If the ing methods generate independently and identically dis- change of measure is well-characterized (as in impor- tributed (i.i.d.) samples. Our method, in contrast, tance sampling), these distortions can be corrected generates independent but non-identically distributed when computing properties of the distribution. (non-i.i.d.) samples. It does so to expose rare behaviors. Our method uses a combination of biased sampling Second, whereas most existing methods for statistical and sequential hypothesis testing [9,10] to estimate the model checking use classical statistics, our method probability that the model satisfies the property. Briefly, employs Bayesian statistics [23,24]. We have previously the algorithm randomly perturbs the computational shown that Bayesian statistics confers advantages in the model prior to generating each sample in order to expose context of statistical verification in terms of efficiency rare but interesting behaviors. These perturbations cause and flexibility [7,8,19,25]. a change of measure. We note that in the context of SDEs, a change of measure is itself a stochastic process Methods and so importance sampling, which assumes that the Our method draws on concepts from several different change of measure with respect to the biased distribution fields. We begin by briefly surveying the semantics of is known, cannot be used. Our technique does not stochastic differential equations, a language for formally require us to know the exact magnitudes of the changes specifying dynamic behaviors, Girsanov’s theorem on of measure, nor the Radon-Nikodym derivatives [11]. change of measures, and results on consistency and con- Instead, it ensures that that the geometric average of centration of Bayesian posteriors. these derivatives is bounded. This is a much weaker assumption than is required by importance sampling, but Stochastic differential equation models we will show that it is sufficient for the purposes of A stochastic differential equation (SDE) [26,27] is a dif- obtaining bounds via sequential hypothesis testing. ferential equation in which some of the terms evolve according to Brownian Motion [28]. A typical SDE is of Related work the following form: Our method performs statistical model checking using dX = b(t, X )dt + v(t, X )dW t t t hypothesis testing [12,13], which has been used pre- viously to analyze in a variety of domains (e.g., where X is a system variable, b is a Riemann integr- [14-19,19]), including computational biology (e.g., able function, v is an Itō integrable function, and W is Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 3 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 Brownian Motion. The Brownian Motion W is a contin- Definition 1 (Adapted Finitely Monitorable). Let s be uous-time stochastic process satisfying the following a finite-length trace from the stochastic differential three conditions: equation . A specification j is said to be adapted finitely monitorable (AFM) if it is possible to decide 1. W =0 whether s satisfies j, denoted j ⊨ j. 2. W is continuous (almost surely). Certain AFM specifications can be expressed as for- 3. W has independent normally distributed mulas in Bounded Linear Temporal Logic (BLTL) increments: [30-32]. Informally, BLTL formulas can capture the � W - W and W - W are independent if 0 ≤ s ordering of events. t s t’ s’ <t <s’ <t’. Definition 2 (Probabilistic Adapted Finitely Monitor- � W − W ∼ N (0, t − s),where N (0, t − s) able). A specification j is said to be probabilistic t s denotes the normal distribution with mean 0 and adapted finitely monitorable (PAFM) if it is possible to variance t - s. Note that the symbol ~ is used to (deterministically or probabilistically) decide whether indicate “is distributed as“. satisfies j with probability at least θ, denoted M  P (φ). ≥θ Consider the time between 0 and t as divided into m Some common examples of PAFM specifications discrete steps 0, t , t ... t = t. The solution to a sto- include Probabilistic Bounded Linear Temporal Logic 1 2 m chastic differential equation is the limit of the following (PBLTL) (e.g., see [19]) and Continuous Specification discrete difference equation, as m goes to infinity: Logic. Note that temporal logic is only one means for constructing specifications; other formalisms can also be X − X = b(t , X )(t − t ) t t k t k+1 k k+1 k k used, like Statecharts [33,34]. + v(t , X )(W − W ) k t t t Semantics of bounded linear temporal logic (BLTL) k k+1 k We define the semantics of BLTL with respect to the In what follows, will refer to a system of stochas- paths of .Let s =(s , Δ ), (s , Δ ),... be a sampled M 0 0 1 1 ticdifferentialequations.Wenotethatasystem of sto- execution of the model along states s , s ,... with dura- 0 1 chastic differential equations comes equipped with an tions Δ , Δ , .. . Î ℝ. We denote the path starting at 0 1 inherent probability space and a natural probability i 0 state i by s (in particular, s denotes the original measure μ. Our algorithm repeatedly and randomly per- execution s). The value of the state variable x in s at turbs the probability measure of the Brownian motion the state i is denoted by V (s, i, x). The semantics of in the model which, in turn, changes the underlying BLTL is defined as follows: measure in an effort to expose rare behaviors. These changes can be characterized using Girsanov’s Theorem. k 1. s ⊨ x ~ v if and only if V(s, k, x)~ v,where v Î Girsanov’s theorem for perturbing stochastic differential ℝ and ~ Î {>, <, =}. equation models k k k 2. s ⊨ j ∨j if and only if s ⊨ j or s ⊨ j . 1 2 1 2 Given a process {θ |0 ≤ t ≤ T}satisfyingthe Novikov k k k 3. s ⊨ j ∧ j if and only if s ⊨ j and s ⊨ j ; 1 2 1 2 condition [29], such as an SDE, the following exponen- k k 4. s ⊨ ¬j if and only if s ⊨ j does not hold. 1 1 tial martingale Z defines the change from measure P to k t 5. s ⊨ j U j if and only if there exists i Î N such 1 2 new measure ˆ : k+i that: (a) 0 ≤ Σ Δ ≤ t;(b) s ⊨ j ;and(c) for 0 ≤ l <i k+l 2 k+j t t each 0 ≤ j<i, s ⊨ j ; Z = exp(− θ dW − θ du/2) t u u 0 0 Statistical model validation Our algorithm performs statistical model checking using Here, Z is the Radon-Nikodym derivative of ˆ with t P Bayesian sequential hypothesis testing [23] on non-i.i.d. respect to P for t <T. The Brownian motion ˆ under ˆ samples. The statistical model checking problem is to ˆ is given by: . The non-stochastic W = W + θ du P t t u decide whether a model satisfies a probabilistic component of the stochastic differential equation is not M adapted finitely monitorable formula j with probability affected by change of measures. Thus, a change of mea- at least θ.Thatis,whether M  P (φ) where θ Î (0, sure for SDEs is a stochastic process (unlike importance ≥θ 1). sampling for explicit probability distributions). Sequential hypothesis testing Let r be the (unknown but fixed) probability of the Specifying dynamic behaviors model satisfying j. We can now re-state the statistical Next, we define a formalism for encoding high-level model checking problem as deciding between the two behavioral specification that our algorithm will test composite hypotheses: H : r ≽ θ and H : r < θ. Here, against . 0 1 M Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 4 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 the null hypothesis H indicates that satisfies the measure of confidence in H vs. H provided by the data 0 0 1 AFM formula j with probability at least θ, while the x ,..., x , and “weighted” by the prior g. 1 n alternate hypothesis H denotes that satisfies the Non-i.i.d. Bayesian sequential hypothesis testing 1 M AFM formula j with probability less than θ. Traditional methods for hypothesis testing, including Definition 3 (Type I and II errors). A Type I error is those outlined in the previous two subsection*s, assume an error where the hypothesis test asserts that the null that the samples are drawn i.i.d.. In this section* we show hypothesis H is false, when in fact H is true. Conver- that non-i.i.d. samples can also be used, provided that cer- 0 0 sely, a Type II error is an error where the hypothesis tain conditions hold. In particular, if one can bound the test asserts that the null hypothesis H is true, when in change in measure associated with the non-identical sam- fact H is false. pling, one can can also bound the Type I and Type II The basic idea behind any statistical model checking errors under a change of measure. Our algorithm bounds algorithm based on sequential hypothesis testing is to the change of measure, and thus the error. iteratively sample traces from the stochastic process. We begin by reviewing some fundamental concepts Each trace is then evaluated with a trace verifier [32], from Bayesian statistics including KL divergence, KL which determines whether the trace satisfies the specifi- support, affinity, and δ-separation, and then restate an cation i.e. s ⊨ j. This is always feasible because the spe- important result on the concentration of Bayesian pos- cifications used are adapted and finitely monitorable. teriors [35,37]. Two accumulators keep track of the total number of Definition 5 (Kullback-Leibler (KL) Divergence). Given traces sampled, and the number of satisfying traces, a parameterized family of probability distributions {f }, respectively. The procedure continues until there is the Kullback-Leibler (KL)divergence K(θ , θ) between enough information to reject either H or H . the distributions corresponding to two parameters θ and 0 1 K (θ , θ)= E [f /f ] Bayesian sequential hypothesis testing θ is: 0 θ θ θ . Note that E is the expecta- 0 θ 0 0 Recall that for any finite trace s of the system and an tion computed under the probability measure f . i 0 Adapted Finitely Monitorable (AFM) formula j,we can Definition 6 (KL Neighborhood). Given a parameter- decide whether s satisfies j. Therefore, we can define a ized family of probability distributions {f }, the KL i θ random variable X denoting the outcome of s ⊨ j. neighborhood K (θ ) of a parameter value θ is given by i i ε 0 0 Thus, X is a Bernoulli random variable with probability the set {θ : K (θ , θ) < ε}. i 0 x 1−x i i mass function ,where x =1 if Definition 7 (KL Support). A point θ is said to be in f (x |ρ)= ρ (1 − ρ) i 0 and only if s ⊨ j, otherwise x =0. the KL support of a prior Π if and only if for all ε >0, Π i q i Bayesian statistics requires that prior probability distri- (K (θ )) >0. ε 0 butions be specified for the unknown quantity, which is Definition 8 (Affinity). The affinity Aff(f, g) between r in our case. Thus we will model the unknown quan- any two densities is defined as . Aff(f , g)= fgdμ tity as a random variable u with prior density g(u). The Definition 9 (Strong δ-Separation). Let A ⊂ [0, 1] and prior probability distribution is usually based on our δ >0. The set A and the point θ are said to be strongly previous experiences and beliefs about the system. Non- δ-separated if and only if for any proper probability dis- informative or objective prior probability distribution tribution v on A, Aff(f , f (x )v(dθ)) <δ . θ θ 1 [35] can be used when nothing is known about the Given these definitions, it can be shown that the Baye- probability of the system satisfying the AFM formula. sian posterior concentrates exponentially under certain Supposewehavea sequence of independent random technical conditions [35,37]. variables X ,..., X defined as above, and let d =(x ,..., Bounding errors under a change of measure 1 n 1 x ) denote a sample of those variables. Next, we develop the machinery needed to compute Definition 4. The Bayes factor of sample d and bounds on the Type-I/Type-II errors for a testing strat- egy based on non-i.i.d. samples. P(d|H ) B = hypotheses H and H is . 0 1 A stochastic differential equation model is natu- P(d|H ) rally associated with a probability measure μ. Our non-i. The Bayes factor may be used as a measure of relative i.d. sampling strategy can be thought of as the assign- confidence in H vs. H , as proposed by Jeffreys [36]. 0 1 ment of a set of probability measures μ , μ ,... to . 1 2 M The Bayes factor is the ratio of two likelihoods: Each unique sample s is associated with an implied 1 probability measure μ and is generated from under i M f (x |u) ··· f (x |u) · g(u)du 1 n μ in an i.i.d. manner. Our proofs require that all the B = . (1) i f (x |u) ··· f (x |u) · g(u)du 1 n implied probability measures are equivalent.Thatis,an event is possible (resp. impossible) under a probability We note that the Bayes factor depends on both the measure if and only if it is possible (resp. impossible) data d and on the prior g, so it may be considered a under the original probability measure μ. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 5 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 We use the following result regarding change of mea- Our result will exploit the fact that we do not allow sures. Suppose a given behavior, say j, holds on the ori- our testing or sampling procedures to have arbitrary ginal model with an (unknown) probability r. implied Radon-Nikodym derivatives. This is reasonable as no statistical guarantees should be available for an p (X |u) g(u) du μ i intelligently designed but adversarial test procedure that P(ρ< θ |X )= 0 i (say) tries to avoid sampling from the given behavior. p (X |u) g(u) du μ i Suppose that the implied Radon-Nikodym derivative Here, X is a Bernoulli random variable denoting the always lies between a constant c and another constant th event that i sample satisfies the given behavior j. Note 1/c. That is, the change of measure does not distort the probabilities of observable events by more than a factor that the X s must be independent of one another. Now, of c. Then, we observe that: we can rewrite the above expression as: p (X |u) θ μ i 0 cp (X |u) g(u) du μ i p (X |u)g(u) du μ i i Q (ρ< θ |X ) ≤ 0 i p (X |u) μ i 1 P (ρ< θ |X )= 0 i p (X |u) g(u) du μ i p (X |u) 1 μ i p (X |u)g(u) du μ i 2 = c P(ρ< θ |X ) 0 i p (X |u) μ i Furthermore, Note that the term p (X |u) denotes the probability of μ i observing the event X under the modified probability Q (ρ< θ |X ) ≥ P (ρ< θ |X ). measure μ if the unknown probability r were u.In 0 i 0 i order to ensure the independence assumption, the new probability measures μ are chosen independently of one Thus, by allowing the sampling algorithm to change p (X |u) μ i another. The ratio is the implied Radon-Niko- measures by at most c,wehavechanged the posterior p (X |u) μ i 2 dym derivative for the change of measure between two probability of observing a behavior by at most c . equivalent probability measures. Suppose, the testing Example: Suppose, the testing strategy has made n strategy has made n observations X , X ,... X . Then, observations X , X ,... X . Then, 1 2 n 1 2 n n n p (X |u)  p (X |u) μ i μ i θ θ i 0 0 p (X |u) g(u) du p (X |u) g(u) du μ i μ i 0 0 p (X |u) p (X |u) μ i μ i i=1 i i=1 P (ρ< θ |X)= Q (ρ< θ |X)= 0 0 n n p (X |u)  p (X |u) μ i μ i 1 1 p (X |u) g(u) du p (X |u) g(u) du μ i μ i 0 0 p (X |u) p (X |u) μ i μ i i=1 i i=1 0 n c (p (X |u)) g(u) du μ i p (X |u) A sampling algorithm can compute by 0 μ i i=1 drawing independent samples from a stochastic differen- tial equation model under the new “modified” probabil- (p (X |u)) g(u) du μ i i=1 ity measure. We note that it is not easy to compute the 2n p (X |u) = c P (ρ< θ |X , X , ... X ) μ i 0 1 2 n change of measure algebraically or numerically. p (X |u) μ i However, our algorithm does not need to compute this Similarly, quantity explicitly. It simply establishes bounds on it. Consider the following expression that is computable (p (X |u)) g(u) du μ i without knowing the implied Radon-Nikodym derivative i=1 Q (ρ< θ |X) ≥ 2n or change of measure explicitly. c (p (X |u)) g(u) du μ i i=1 p (X |u) g(u) du μ i 0 i Q (ρ< θ |X )= = P (ρ< θ |X , ... X ) 0 i  0 1 n 2n p (X |u) g(u) du μ i Termination conditions for non-i.i.d. sampling Now, we can rewrite the above expression as: Traditional (i.e., i.i.d.) Bayesian Sequential Hypothesis p (X |u) Testing is guaranteed to terminate. That is, only a finite θ μ i 0 i p (X |u) g(u) du μ i number of samples are required before the test selects p (X |u) μ i Q (ρ< θ |X )= 0 i one of the hypotheses. We now consider the conditions p (X |u) μ i p (X |u) g(u) du μ i 0 under which a Bayesian Sequential Hypothesis Testing p (X |u) μ i based procedure using non-i.i.d. samples will terminate. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 6 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 To do this, we first need to show that the posterior denoted the number of samples that satisfy the AFM probability distribution will concentrate on a particular specification j. Based upon the samples observed, we value as we see more an more samples from the model. compute the Bayes Factor under the new probability To consider the conditions under which our algorithm measures. We know that the Bayes Factor so computed will terminate after observing n samples, note that the is within a factor of the original Bayes Factor under 2n factor introduced due to the change of measure c can the natural probability measure. Hence, the algorithm 2n outweigh the gain made by the concentration of the divides the Bayes Factor by the factor h if the Bayes -nb probability measure e . This is not surprising because Factor is larger than one. If the Bayes Factor is less our construction thus far does not force the test not to than one, the algorithm multiplies the Bayes Factor by 2n bias against a sample in an intelligent way. That is, a the factor h . maliciously designed testing procedure could simply avoid the error prone regions of the design. To address Results and discussion this, we define the notion of a fair testing strategy that We applied our algorithm to two SDE models of tumor does not engage in such malicious sampling. dynamics from the literature. The first model is a single Definition 10. A testing strategy is h-fair (h ≥ 1) if dimensional stochastic differential equation for the and only if the geometric average of the implied Radon- influence of chemotherapy on cancer cells, and the sec- Nikodym derivatives over a number of samples is within ond model is a pair of SDEs that describe an immuno- a constant factor h of unity, i.e., genic tumor. Lefever and Garay model 1 p (X |u) μ i ≤ ≤ η Lefever and Garay [38] studied the growth of tumors η p (X |u) μ i i=1 under immune surveillance and chemotherapy using the following stochastic differential equation: Note that a fair test strategy does not need to sample from the underlying distribution in an i.i.d. manner. dx x βx x = r x (1 − ) − + x(1 − ) A cos (ωt) 0 0 However, it must guarantee that the probability of dt K 1+ x K observing the given behavior in a large number of + x (1 − ) W observations is not altered substantially by the non-i.i.d. sampling. Intuitively, we want to make sure that we bias Here, x is the amount of tumor cells, A cos(ωt) for each sample as many times as we bias against it. denotes the influence of a periodic chemotherapy treat- Our main result shows that such a long term neutrality ment, r is the linear per capita birth rate of cancer is sufficient to generate statistical guarantees on an cells, K is the carrying capacity of the environment, and otherwise non-i.i.d. testing procedure. b represents the influence of the immune cells. Note Definition 11. An h-fair test is said to be eventually 4 b that W is the standard Brownian Motion, and default fair if and only if 1 ≤ h <e ,where b is the constant in model parameters were those used in [38].. the exponential posterior concentration theorem. We demonstrate our algorithm on a simple property The notion of a eventually fair test corresponds to a of the model. Namely, starting with a tumor consisting testing strategy that is not malicious or adversarial, and of a billion cells, is there at least a 1% chance that the is making an honest attempt to sample from all the tumor could increase to one hundred billion cells under events in the long run. under immune surveillance and chemotherapy. The fol- lowing BLTL specification captures the behavioral speci- Algorithm fication: Finally, we present our Statistical Verification algo- rithm(SeeFigure2)interms of ageneric non-i.i.d. 10 11 Pr (F (x> 10 )) ≥0.01 testing procedure sampling with random “implied” change of measures. Our algorithm is relatively simple Figure 3 contrasts the number of samples needed to and generalizes our previous Bayesian Statistical verifi- decide whether the model satisfies the property using i.i. cation algorithm [8] to non-i.i.d. samples using change d. and non-i.i.d.. As expected, the number of samples of measures. The algorithm draws non-i.i.d. samples required increases linearly in the logarithm of the Bayes from the stochastic differential equation under ran- factor regardless of whether i.i.d.ornon-i.i.d. sampling domly chosen probability measures. The algorithm is used. However, non-i.i.d. sampling always requires ensures that the implied change of measure is bounded fewer samples than i.i.d. sampling. Moreover, the differ- so as to make the testing approach fair. The variable n ence between the number of samples increases with the denotes the number of samples obtained so far and x Bayes factor. That is, the lines are diverging. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 7 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 Figure 2 Non-i.i.d. Statistical Verification Algorithm. The figure illustrated the non-i.i.d. Bayesian model validation algorithm. The algorithm builds upon Girsanov’s theorem on change of measure and Bayesian model validation. Figure 3 Comparison of i.i.d. and non-i.i.d. sampling. Non-i.i.d. vs i.i.d. sampling based verification for the Lefever and Garay model. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 8 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 We note that there are circumstances when our algo- For this model, we considered the following property: rithm may require more samples than one based on i.i. starting from 0.1 units each of tumor and immune cells, d. sampling. This will happen when the property being is thereatleast a1%chancethatthe number of tumor tested has relatively high probability. For example, we cells could increase to 3.3 units. The property can be tested the property that the probability of eradicating encoded into the following BLTL specification: the tumor is at most 1%. The i.i.d. algorithm required 1374 samples to deny this possibility while the non-i.i.d. Pr (F (x > 3.3)) ≥0.01 algorithm required 1526 samples. Thus, our algorithm is best used to examine rare behaviors. Default model parameters were those used in [39]. Fig- ure 4 contrasts the number of samples needed to decide Nonlinear immunogenic tumor model whether the model satisfies the property using i.i.d.and The second model we analyze studies immunogenic non-i.i.d.. The same trends are observed as in the pre- tumor growth [39,40]. Unlike the previous model, the vious model. That is, our algorithm requires fewer sam- immunogenic tumor model explicitly tracks the ples than i.i.d. hypothesis testing, and that the difference dynamics of the immune cells (variable x) in response to between these methods increases with the Bayes factor. the tumor cells (variable y). The SDEs are as follows: We also considered the property that the number of tumor cells increases to 4.0 units. We evaluated whether dx(t)= (a − a x (t)+ a x (t)y(t)) dt 1 2 3 this property is true with probability at least 0.000005 +(b (x(t) − x )+ b (y(t) − y )) dW (t) 11 1 12 1 1 under a Bayes Factor of 100, 000. The i.i.d.sampling dy(t)= (b y(t)(1 − b y(t)) − x(t)y(t)) dt 1 2 algorithm did not produce an answer even after obser- +(b (x(t) − x )+ b (y(t) − y )) dW (t) 21 1 22 1 2 ving 10, 000 samples. The non-i.i.d. model validation algorithm answered affirmatively after observing 6, 736 The parameters x an y denote the stochastic equili- 1 1 samples. Once again, the real impact of the proposed brium point of the model. Briefly, the model assumes algorithm lies in uncovering rare behaviors and bound- that the amount of noise increases with the distance to ing their probability of occurrence. the equilibrium point. Figure 4 Comparison of i.i.d. and non-i.i.d. sampling. Non-i.i.d. vs i.i.d. Sampling based verification for the nonlinear Immunogenic tumor model. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 9 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 References Discussion 1. Faeder JR, Blinov ML, Goldstein B, Hlavacek WS: Rule-based modeling of Our results confirm that non-i.i.d. sampling reduces the biochemical networks. Complexity 2005, 10(4):22-41. number of samples required in the context of hypothesis 2. Haigh J: Stochastic Modelling for Systems Biology by D. J. Wilkinson. Journal Of The Royal Statistical Society Series A 2007, 170:261-261[http:// testing – when the property under consideration is rare. ideas.repec. org/a/bla/jorssa/v170y2007i1p261-261.html]. Moreover, the benefits of non-i.i.d. sampling increase 3. Iosifescu M, Tautu P, Iosifescu M: Stochastic processes and applications in with the rarity of the property, as confirmed by the biology and medicine [by] M. Iosifescu [and] P. Tautu Editura Academiei; Springer-Verlag, Bucuresti, New York; 1973. divergence of the lines in Figures 2 and 3. We note, 4. Twycross J, Band L, Bennett MJ, King J, Krasnogor N: Stochastic and however, that if the property isn’trarethenanon-i.i.d. deterministic multiscale models for systems biology: an auxin-transport sampling strategy will actually require a larger number case study. BMC Systems Biology 2010, 4(34)[http://www. biomedcentral. com/17520509/4/34/abstract]. of samples than an i.i.d. strategy. Thus, our algorithm is 5. Kwiatkowska M, Norman G, Parker D: Stochastic Model Checking. In only appropriate for investigating rare behaviors. Formal Methods for the Design of Computer, Communication and Software Systems: Performance Evaluation (SFM’07). Volume 4486. Springer;Bernardo M, Hillston J 2007:220-270, LNCS (Tutorial Volume). Conclusions 6. Kwiatkowska M, Norman G, Parker D: Advances and Challenges of We have introduced the first algorithm for verifying Probabilistic Model Checking. Proc 48th Annual Allerton Conference on properties of stochastic differential equations using Communication, Control and Computing IEEE Press; 2010. 7. Langmead C: Generalized Queries and Bayesian Statistical Model sequential hypothesis testing. Our technique combines Checking in Dynamic Bayesian Networks: Application to Personalized Bayesian statistical model checking and non-i.i.d.sam- Medicine. Proc of the 8th International Conference on Computational Systems pling and provides guarantees in terms of termination, Bioinformatics (CSB) 2009, 201-212. 8. Jha SK, Clarke EM, Langmead CJ, Legay A, Platzer A, Zuliani P: A Bayesian and the number of samples needed to achieve those Approach to Model Checking Biological Systems. In CMSB Volume 5688 of bounds. The method is most suitable when the behavior Lecture Notes in Computer Science. Springer;Degano P, Gorrieri R of interest is the exception and not the norm. 2009:218-234. 9. Wald A: Sequential Analysis New York: John Wiley and Son; 1947. The present paper only considers SDEs with indepen- 10. Lehmann EL, Romano JP: Testing statistical hypotheses. 3 edition. Springer dent Brownian noise. We believe that these results can Texts in Statistics, New York: Springer; 2005. be extended to handle SDEs with certain kinds of corre- 11. Edgar GA: Radon-Nikodym theorem. Duke Mathematical Journal 1975, 42:447-450. lated noise. Another interesting direction for future 12. Younes HLS, Simmons RG: Probabilistic Verification of Discrete Event work is the extension of these method to stochastic par- Systems Using Acceptance Sampling. In CAV, Volume 2404 of Lecture Notes tial differential equations, which are used to model spa- in Computer Science. Springer;Brinksma E, Larsen KG 2002:223-235. 13. Younes HLS, Kwiatkowska MZ, Norman G, Parker D: Numerical vs. tially inhomogeneous processes. Such analysis methods Statistical Probabilistic Model Checking: An Empirical Study. TACAS 2004, could be used, for example, to investigate properties 46-60. concerning spatial properties of tumors, the propagation 14. Lassaigne R, Peyronnet S: Approximate Verification of Probabilistic Systems. PAPM-PROBMIV 2002, 213-214. of electrical waves in cardiac tissue, or more generally, 15. Hérault T, Lassaigne R, Magniette F, Peyronnet S: Approximate Probabilistic to the diffusion processes observed in nature. Model Checking. Proc 5th International Conference on Verification, Model Checking and Abstract Interpretation (VMCAI’04), Volume 2937 of LNCS Springer; 2004. Acknowledgements 16. Sen K, Viswanathan M, Agha G: Statistical Model Checking of Black-Box The authors acknowledge the feedback received from the anonymous Probabilistic Systems. CAV 2004, 202-215. reviewers for the first IEEE Conference on Compuational Advances in Bio 17. Grosu R, Smolka S: Monte Carlo Model Checking. CAV 2005, 271-286. and Medical Sciences (ICCABS) 2011. 18. Gondi K, Patel Y, Sistla AP: Monitoring the Full Range of omega-Regular This article has been published as part of BMC Bioinformatics Volume 13 Properties of Stochastic Systems. Verification, Model Checking, and Abstract Supplement 5, 2012: Selected articles from the First IEEE International Interpretation, 10th International Conference, VMCAI 2009, Volume 5403 of Conference on Computational Advances in Bio and medical Sciences LNCS Springer; 2009, 105-119. (ICCABS 2011): Bioinformatics. The full contents of the supplement are 19. Jha SK, Langmead C, Ramesh S, Mohalik S: When to stop verification? available online at http://www.biomedcentral.com/bmcbioinformatics/ Statistical Trade-off between Costs and Expected Losses. Proceedings of supplements/13/S5. Design Automation and test In Europe (DATE) 2011, 1309-1314. 20. Clarke EM, Faeder JR, Langmead CJ, Harris LA, Jha SK, Legay A: Statistical Author details Model Checking in BioLab: Applications to the Automated Analysis of T- Electrical Engineering and Computer Science Department, University of Cell Receptor Signaling Pathway. CMSB 2008, 231-250. Central Florida, Orlando FL 32816 USA. Computer Science Department, 21. Langmead CJ, Jha SK: Predicting Protein Folding Kinetics Via Temporal Carnegie Mellon University, Pittsburgh PA 15213 USA. Lane Center for Logic Model Checking. In WABI, Volume 4645 of Lecture Notes in Computer Computational Biology, Carnegie Mellon University, Pittsburgh PA 15213 Science. Springer;Giancarlo R, Hannenhalli S 2007:252-264. USA. 22. Jha S, Langmead C: Synthesis and Infeasibility Analysis for Stochastic Models of Biochemical Systems using Statistical Model Checking and Authors’ contributions Abstraction Refinement.. SKJ and CJL contributed equally to all parts of the paper. 23. Jeffreys H: Theory of Probability Oxford University Press; 1939. 24. Gelman A, Carlin JB, Stern HS, Rubin DB: Bayesian Data Analysis London: Competing interests Chapman & Hall; 1995. The authors declare that they have no competing interests. 25. Jha S, Langmead C: Exploring Behaviors of SDE Models of Biological Systems using Change of Measures. Proc of the 1st IEEE International Published: 12 April 2012 Conference on Computational Advances in Bio and medical Sciences (ICCABS) 2011, 111-117. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 10 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 26. Øksendal B: Stochastic Differential Equations: An Introduction with Applications (Universitext). 6 edition. Springer; 2003 [http://www.amazon. com/exec/obidos/redirect?tag=citeulike07-20\&path=ASIN/3540047581]. 27. Karatzas I, Shreve SE: Brownian Motion and Stochastic Calculus (Graduate Texts in Mathematics). 2 edition. Springer; 1991 [http://www.amazon.com/ exec/obidos/redirect?tag=citeulike07-20\&path=ASIN/0387976558]. 28. Wang MC, Uhlenbeck GE: On the Theory of the Brownian Motion II. Reviews of Modern Physics 1945, 17(2-3):323[http://dx.doi.org/10.1103/ RevModPhys. 17.323]. 29. Girsanov IV: On Transforming a Certain Class of Stochastic Processes by Absolutely Continuous Substitution of Measures. Theory of Probability and its Applications 1960, 5(3):285-301. 30. Pnueli A: The Temporal Logic of Programs. FOCS IEEE; 1977, 46-57. 31. Owicki SS, Lamport L: Proving Liveness Properties of Concurrent Programs. ACM Trans Program Lang Syst 1982, 4(3):455-495. 32. Finkbeiner B, Sipma H: Checking Finite Traces Using Alternating Automata. Formal Methods in System Design 2004, 24(2):101-127. 33. Harel D: Statecharts: A Visual Formalism for Complex Systems. Science of Computer Programming 1987, 8(3):231-274[http://citeseerx.ist.psu.edu/ viewdoc/summary?doi=10.1.1.20.4799]. 34. Ehrig H, Orejas F, Padberg J: Relevance, Integration and Classification of Specification Formalisms and Formal Specification Techniques.[http:// citeseerx.ist. psu.edu/viewdoc/summary?doi=10.1.1.42.7137]. 35. Berger J: Statistical Decision Theory and Bayesian Analysis Springer-Verlag; 36. Jeffreys H: Theory of probability/by Harold Jeffreys. 3 edition. Clarendon Press, Oxford; 1961. 37. Choi T, Ramamoorthi RV: Remarks on consistency of posterior distributions. ArXiv e-prints 2008. 38. Lefever R, RGaray S: Biomathematics and Cell Kinetics Elsevier, North-Hollan biomedical Press; 1978, chap. Local description of immune tumor rejection,:333. 39. Horhat R, Horhat R, Opris D: The simulation of a stochastic model for tumour-immune system. Proceedings of the 2nd WSEAS international conference on Biomedical electronics and biomedical informatics BEBI’09, Stevens Point, Wisconsin, USA: World Scientific and Engineering Academy and Society (WSEAS); 2009, 247-252[http://portal.acm.org/citation.cfm?id = 1946539.1946584]. 40. Kuznetsov V, Makalkin I, Taylor M, Perelson A: Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bulletin of Mathematical Biology 1994, 56:295-321[http://dx. doi. org/10.1007/BF02460644], 10.1007/BF02460644. doi:10.1186/1471-2105-13-S5-S8 Cite this article as: Jha and Langmead: Exploring behaviors of stochastic differential equation models of biological systems using change of measures. BMC Bioinformatics 2012 13(Suppl 5):S8. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png BMC Bioinformatics Springer Journals

Exploring behaviors of stochastic differential equation models of biological systems using change of measures

Jha, Sumit; Langmead, Christopher
BMC Bioinformatics , Volume 13 (5) – Apr 12, 2012
Download PDF
10 pages

Loading next page...
 
/lp/springer-journals/exploring-behaviors-of-stochastic-differential-equation-models-of-GlaG21tqM0

References (76)

Publisher
Springer Journals
Copyright
Copyright © 2012 by Jha and Langmead; licensee BioMed Central Ltd.
Subject
Life Sciences; Bioinformatics; Microarrays; Computational Biology/Bioinformatics; Computer Appl. in Life Sciences; Combinatorial Libraries; Algorithms
eISSN
1471-2105
DOI
10.1186/1471-2105-13-S5-S8
pmid
22537012
Publisher site
See Article on Publisher Site

Abstract

Stochastic Differential Equations (SDE) are often used to model the stochastic dynamics of biological systems. Unfortunately, rare but biologically interesting behaviors (e.g., oncogenesis) can be difficult to observe in stochastic models. Consequently, the analysis of behaviors of SDE models using numerical simulations can be challenging. We introduce a method for solving the following problem: given a SDE model and a high-level behavioral specification about the dynamics of the model, algorithmically decide whether the model satisfies the specification. While there are a number of techniques for addressing this problem for discrete-state stochastic models, the analysis of SDE and other continuous-state models has received less attention. Our proposed solution uses a combination of Bayesian sequential hypothesis testing, non-identically distributed samples, and Girsanov’s theorem for change of measures to examine rare behaviors. We use our algorithm to analyze two SDE models of tumor dynamics. Our use of non-identically distributed samples sampling contributes to the state of the art in statistical verification and model checking of stochastic models by providing an effective means for exposing rare events in SDEs, while retaining the ability to compute bounds on the probability that those events occur. Background level description of a dynamical behavior, j (e.g., a for- The dynamics of biological systems are largely driven by mula in temporal logic). It then computes a statistically stochastic processes and subject to random external per- rigorous bound on the probability that the given model turbations. The consequences of such random processes exhibits the stated behavior using a combination of are often investigated through the development and ana- biased sampling and Bayesian Statistical Model Check- lysis of stochastic models (e.g., [1-4]). Unfortunately, the ing [7,8]. validation and analysis of stochastic models can be very Existing methods for validating and analyzing stochas- challenging [5,6], especially when the model is intended tic models often require extensive Monte Carlo sam- to investigate rare, but biologically significant behaviors pling of independent trajectories to verify that the (e.g., oncogenesis). The goal of this paper is to introduce model is consistent with known data, and to character- an algorithm for examining such rare behaviors in Sto- ize the model’s expected behavior under various initial chastic Differential Equation (SDE) models. The algo- conditions. Sampling strategies are either unbiased or rithm takes as input the SDE model, ,and ahigh- biased. Unbiased sampling strategies draw trajectories according to the probability distribution implied by the model, and are thus not well-suited to investigating rare * Correspondence: [email protected]; [email protected] 1 behaviors. For example, if the actual probability that the Electrical Engineering and Computer Science Department, University of -10 Central Florida, Orlando FL 32816 USA model will exhibit a given behavior is 10 , then the Computer Science Department, Carnegie Mellon University, Pittsburgh PA expected number of samples need to see such behaviors 15213 USA is about 10 (See Figure 1). Biased sampling strategies, Full list of author information is available at the end of the article © 2012 Jha and Langmead; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 2 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 Figure 1 Observing rare behaviors in i.i.d. sampling is challenging. A toy model with a one low-probability state. An unbiased sampling algorithm may require billions of samples in order to observe the ‘bad’ state. Statistical algorithms based on i.i.d. algorithms are not suitable for analyzing such models with rare interesting behaviors. in contrast, can be used to increase the probability of [7,8,20-22]). There are two key differences between the observing rare events, at the expense of distorting the method in this paper and existing methods. First, exist- underlying probability distribution/measure. If the ing methods generate independently and identically dis- change of measure is well-characterized (as in impor- tributed (i.i.d.) samples. Our method, in contrast, tance sampling), these distortions can be corrected generates independent but non-identically distributed when computing properties of the distribution. (non-i.i.d.) samples. It does so to expose rare behaviors. Our method uses a combination of biased sampling Second, whereas most existing methods for statistical and sequential hypothesis testing [9,10] to estimate the model checking use classical statistics, our method probability that the model satisfies the property. Briefly, employs Bayesian statistics [23,24]. We have previously the algorithm randomly perturbs the computational shown that Bayesian statistics confers advantages in the model prior to generating each sample in order to expose context of statistical verification in terms of efficiency rare but interesting behaviors. These perturbations cause and flexibility [7,8,19,25]. a change of measure. We note that in the context of SDEs, a change of measure is itself a stochastic process Methods and so importance sampling, which assumes that the Our method draws on concepts from several different change of measure with respect to the biased distribution fields. We begin by briefly surveying the semantics of is known, cannot be used. Our technique does not stochastic differential equations, a language for formally require us to know the exact magnitudes of the changes specifying dynamic behaviors, Girsanov’s theorem on of measure, nor the Radon-Nikodym derivatives [11]. change of measures, and results on consistency and con- Instead, it ensures that that the geometric average of centration of Bayesian posteriors. these derivatives is bounded. This is a much weaker assumption than is required by importance sampling, but Stochastic differential equation models we will show that it is sufficient for the purposes of A stochastic differential equation (SDE) [26,27] is a dif- obtaining bounds via sequential hypothesis testing. ferential equation in which some of the terms evolve according to Brownian Motion [28]. A typical SDE is of Related work the following form: Our method performs statistical model checking using dX = b(t, X )dt + v(t, X )dW t t t hypothesis testing [12,13], which has been used pre- viously to analyze in a variety of domains (e.g., where X is a system variable, b is a Riemann integr- [14-19,19]), including computational biology (e.g., able function, v is an Itō integrable function, and W is Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 3 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 Brownian Motion. The Brownian Motion W is a contin- Definition 1 (Adapted Finitely Monitorable). Let s be uous-time stochastic process satisfying the following a finite-length trace from the stochastic differential three conditions: equation . A specification j is said to be adapted finitely monitorable (AFM) if it is possible to decide 1. W =0 whether s satisfies j, denoted j ⊨ j. 2. W is continuous (almost surely). Certain AFM specifications can be expressed as for- 3. W has independent normally distributed mulas in Bounded Linear Temporal Logic (BLTL) increments: [30-32]. Informally, BLTL formulas can capture the � W - W and W - W are independent if 0 ≤ s ordering of events. t s t’ s’ <t <s’ <t’. Definition 2 (Probabilistic Adapted Finitely Monitor- � W − W ∼ N (0, t − s),where N (0, t − s) able). A specification j is said to be probabilistic t s denotes the normal distribution with mean 0 and adapted finitely monitorable (PAFM) if it is possible to variance t - s. Note that the symbol ~ is used to (deterministically or probabilistically) decide whether indicate “is distributed as“. satisfies j with probability at least θ, denoted M  P (φ). ≥θ Consider the time between 0 and t as divided into m Some common examples of PAFM specifications discrete steps 0, t , t ... t = t. The solution to a sto- include Probabilistic Bounded Linear Temporal Logic 1 2 m chastic differential equation is the limit of the following (PBLTL) (e.g., see [19]) and Continuous Specification discrete difference equation, as m goes to infinity: Logic. Note that temporal logic is only one means for constructing specifications; other formalisms can also be X − X = b(t , X )(t − t ) t t k t k+1 k k+1 k k used, like Statecharts [33,34]. + v(t , X )(W − W ) k t t t Semantics of bounded linear temporal logic (BLTL) k k+1 k We define the semantics of BLTL with respect to the In what follows, will refer to a system of stochas- paths of .Let s =(s , Δ ), (s , Δ ),... be a sampled M 0 0 1 1 ticdifferentialequations.Wenotethatasystem of sto- execution of the model along states s , s ,... with dura- 0 1 chastic differential equations comes equipped with an tions Δ , Δ , .. . Î ℝ. We denote the path starting at 0 1 inherent probability space and a natural probability i 0 state i by s (in particular, s denotes the original measure μ. Our algorithm repeatedly and randomly per- execution s). The value of the state variable x in s at turbs the probability measure of the Brownian motion the state i is denoted by V (s, i, x). The semantics of in the model which, in turn, changes the underlying BLTL is defined as follows: measure in an effort to expose rare behaviors. These changes can be characterized using Girsanov’s Theorem. k 1. s ⊨ x ~ v if and only if V(s, k, x)~ v,where v Î Girsanov’s theorem for perturbing stochastic differential ℝ and ~ Î {>, <, =}. equation models k k k 2. s ⊨ j ∨j if and only if s ⊨ j or s ⊨ j . 1 2 1 2 Given a process {θ |0 ≤ t ≤ T}satisfyingthe Novikov k k k 3. s ⊨ j ∧ j if and only if s ⊨ j and s ⊨ j ; 1 2 1 2 condition [29], such as an SDE, the following exponen- k k 4. s ⊨ ¬j if and only if s ⊨ j does not hold. 1 1 tial martingale Z defines the change from measure P to k t 5. s ⊨ j U j if and only if there exists i Î N such 1 2 new measure ˆ : k+i that: (a) 0 ≤ Σ Δ ≤ t;(b) s ⊨ j ;and(c) for 0 ≤ l <i k+l 2 k+j t t each 0 ≤ j<i, s ⊨ j ; Z = exp(− θ dW − θ du/2) t u u 0 0 Statistical model validation Our algorithm performs statistical model checking using Here, Z is the Radon-Nikodym derivative of ˆ with t P Bayesian sequential hypothesis testing [23] on non-i.i.d. respect to P for t <T. The Brownian motion ˆ under ˆ samples. The statistical model checking problem is to ˆ is given by: . The non-stochastic W = W + θ du P t t u decide whether a model satisfies a probabilistic component of the stochastic differential equation is not M adapted finitely monitorable formula j with probability affected by change of measures. Thus, a change of mea- at least θ.Thatis,whether M  P (φ) where θ Î (0, sure for SDEs is a stochastic process (unlike importance ≥θ 1). sampling for explicit probability distributions). Sequential hypothesis testing Let r be the (unknown but fixed) probability of the Specifying dynamic behaviors model satisfying j. We can now re-state the statistical Next, we define a formalism for encoding high-level model checking problem as deciding between the two behavioral specification that our algorithm will test composite hypotheses: H : r ≽ θ and H : r < θ. Here, against . 0 1 M Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 4 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 the null hypothesis H indicates that satisfies the measure of confidence in H vs. H provided by the data 0 0 1 AFM formula j with probability at least θ, while the x ,..., x , and “weighted” by the prior g. 1 n alternate hypothesis H denotes that satisfies the Non-i.i.d. Bayesian sequential hypothesis testing 1 M AFM formula j with probability less than θ. Traditional methods for hypothesis testing, including Definition 3 (Type I and II errors). A Type I error is those outlined in the previous two subsection*s, assume an error where the hypothesis test asserts that the null that the samples are drawn i.i.d.. In this section* we show hypothesis H is false, when in fact H is true. Conver- that non-i.i.d. samples can also be used, provided that cer- 0 0 sely, a Type II error is an error where the hypothesis tain conditions hold. In particular, if one can bound the test asserts that the null hypothesis H is true, when in change in measure associated with the non-identical sam- fact H is false. pling, one can can also bound the Type I and Type II The basic idea behind any statistical model checking errors under a change of measure. Our algorithm bounds algorithm based on sequential hypothesis testing is to the change of measure, and thus the error. iteratively sample traces from the stochastic process. We begin by reviewing some fundamental concepts Each trace is then evaluated with a trace verifier [32], from Bayesian statistics including KL divergence, KL which determines whether the trace satisfies the specifi- support, affinity, and δ-separation, and then restate an cation i.e. s ⊨ j. This is always feasible because the spe- important result on the concentration of Bayesian pos- cifications used are adapted and finitely monitorable. teriors [35,37]. Two accumulators keep track of the total number of Definition 5 (Kullback-Leibler (KL) Divergence). Given traces sampled, and the number of satisfying traces, a parameterized family of probability distributions {f }, respectively. The procedure continues until there is the Kullback-Leibler (KL)divergence K(θ , θ) between enough information to reject either H or H . the distributions corresponding to two parameters θ and 0 1 K (θ , θ)= E [f /f ] Bayesian sequential hypothesis testing θ is: 0 θ θ θ . Note that E is the expecta- 0 θ 0 0 Recall that for any finite trace s of the system and an tion computed under the probability measure f . i 0 Adapted Finitely Monitorable (AFM) formula j,we can Definition 6 (KL Neighborhood). Given a parameter- decide whether s satisfies j. Therefore, we can define a ized family of probability distributions {f }, the KL i θ random variable X denoting the outcome of s ⊨ j. neighborhood K (θ ) of a parameter value θ is given by i i ε 0 0 Thus, X is a Bernoulli random variable with probability the set {θ : K (θ , θ) < ε}. i 0 x 1−x i i mass function ,where x =1 if Definition 7 (KL Support). A point θ is said to be in f (x |ρ)= ρ (1 − ρ) i 0 and only if s ⊨ j, otherwise x =0. the KL support of a prior Π if and only if for all ε >0, Π i q i Bayesian statistics requires that prior probability distri- (K (θ )) >0. ε 0 butions be specified for the unknown quantity, which is Definition 8 (Affinity). The affinity Aff(f, g) between r in our case. Thus we will model the unknown quan- any two densities is defined as . Aff(f , g)= fgdμ tity as a random variable u with prior density g(u). The Definition 9 (Strong δ-Separation). Let A ⊂ [0, 1] and prior probability distribution is usually based on our δ >0. The set A and the point θ are said to be strongly previous experiences and beliefs about the system. Non- δ-separated if and only if for any proper probability dis- informative or objective prior probability distribution tribution v on A, Aff(f , f (x )v(dθ)) <δ . θ θ 1 [35] can be used when nothing is known about the Given these definitions, it can be shown that the Baye- probability of the system satisfying the AFM formula. sian posterior concentrates exponentially under certain Supposewehavea sequence of independent random technical conditions [35,37]. variables X ,..., X defined as above, and let d =(x ,..., Bounding errors under a change of measure 1 n 1 x ) denote a sample of those variables. Next, we develop the machinery needed to compute Definition 4. The Bayes factor of sample d and bounds on the Type-I/Type-II errors for a testing strat- egy based on non-i.i.d. samples. P(d|H ) B = hypotheses H and H is . 0 1 A stochastic differential equation model is natu- P(d|H ) rally associated with a probability measure μ. Our non-i. The Bayes factor may be used as a measure of relative i.d. sampling strategy can be thought of as the assign- confidence in H vs. H , as proposed by Jeffreys [36]. 0 1 ment of a set of probability measures μ , μ ,... to . 1 2 M The Bayes factor is the ratio of two likelihoods: Each unique sample s is associated with an implied 1 probability measure μ and is generated from under i M f (x |u) ··· f (x |u) · g(u)du 1 n μ in an i.i.d. manner. Our proofs require that all the B = . (1) i f (x |u) ··· f (x |u) · g(u)du 1 n implied probability measures are equivalent.Thatis,an event is possible (resp. impossible) under a probability We note that the Bayes factor depends on both the measure if and only if it is possible (resp. impossible) data d and on the prior g, so it may be considered a under the original probability measure μ. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 5 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 We use the following result regarding change of mea- Our result will exploit the fact that we do not allow sures. Suppose a given behavior, say j, holds on the ori- our testing or sampling procedures to have arbitrary ginal model with an (unknown) probability r. implied Radon-Nikodym derivatives. This is reasonable as no statistical guarantees should be available for an p (X |u) g(u) du μ i intelligently designed but adversarial test procedure that P(ρ< θ |X )= 0 i (say) tries to avoid sampling from the given behavior. p (X |u) g(u) du μ i Suppose that the implied Radon-Nikodym derivative Here, X is a Bernoulli random variable denoting the always lies between a constant c and another constant th event that i sample satisfies the given behavior j. Note 1/c. That is, the change of measure does not distort the probabilities of observable events by more than a factor that the X s must be independent of one another. Now, of c. Then, we observe that: we can rewrite the above expression as: p (X |u) θ μ i 0 cp (X |u) g(u) du μ i p (X |u)g(u) du μ i i Q (ρ< θ |X ) ≤ 0 i p (X |u) μ i 1 P (ρ< θ |X )= 0 i p (X |u) g(u) du μ i p (X |u) 1 μ i p (X |u)g(u) du μ i 2 = c P(ρ< θ |X ) 0 i p (X |u) μ i Furthermore, Note that the term p (X |u) denotes the probability of μ i observing the event X under the modified probability Q (ρ< θ |X ) ≥ P (ρ< θ |X ). measure μ if the unknown probability r were u.In 0 i 0 i order to ensure the independence assumption, the new probability measures μ are chosen independently of one Thus, by allowing the sampling algorithm to change p (X |u) μ i another. The ratio is the implied Radon-Niko- measures by at most c,wehavechanged the posterior p (X |u) μ i 2 dym derivative for the change of measure between two probability of observing a behavior by at most c . equivalent probability measures. Suppose, the testing Example: Suppose, the testing strategy has made n strategy has made n observations X , X ,... X . Then, observations X , X ,... X . Then, 1 2 n 1 2 n n n p (X |u)  p (X |u) μ i μ i θ θ i 0 0 p (X |u) g(u) du p (X |u) g(u) du μ i μ i 0 0 p (X |u) p (X |u) μ i μ i i=1 i i=1 P (ρ< θ |X)= Q (ρ< θ |X)= 0 0 n n p (X |u)  p (X |u) μ i μ i 1 1 p (X |u) g(u) du p (X |u) g(u) du μ i μ i 0 0 p (X |u) p (X |u) μ i μ i i=1 i i=1 0 n c (p (X |u)) g(u) du μ i p (X |u) A sampling algorithm can compute by 0 μ i i=1 drawing independent samples from a stochastic differen- tial equation model under the new “modified” probabil- (p (X |u)) g(u) du μ i i=1 ity measure. We note that it is not easy to compute the 2n p (X |u) = c P (ρ< θ |X , X , ... X ) μ i 0 1 2 n change of measure algebraically or numerically. p (X |u) μ i However, our algorithm does not need to compute this Similarly, quantity explicitly. It simply establishes bounds on it. Consider the following expression that is computable (p (X |u)) g(u) du μ i without knowing the implied Radon-Nikodym derivative i=1 Q (ρ< θ |X) ≥ 2n or change of measure explicitly. c (p (X |u)) g(u) du μ i i=1 p (X |u) g(u) du μ i 0 i Q (ρ< θ |X )= = P (ρ< θ |X , ... X ) 0 i  0 1 n 2n p (X |u) g(u) du μ i Termination conditions for non-i.i.d. sampling Now, we can rewrite the above expression as: Traditional (i.e., i.i.d.) Bayesian Sequential Hypothesis p (X |u) Testing is guaranteed to terminate. That is, only a finite θ μ i 0 i p (X |u) g(u) du μ i number of samples are required before the test selects p (X |u) μ i Q (ρ< θ |X )= 0 i one of the hypotheses. We now consider the conditions p (X |u) μ i p (X |u) g(u) du μ i 0 under which a Bayesian Sequential Hypothesis Testing p (X |u) μ i based procedure using non-i.i.d. samples will terminate. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 6 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 To do this, we first need to show that the posterior denoted the number of samples that satisfy the AFM probability distribution will concentrate on a particular specification j. Based upon the samples observed, we value as we see more an more samples from the model. compute the Bayes Factor under the new probability To consider the conditions under which our algorithm measures. We know that the Bayes Factor so computed will terminate after observing n samples, note that the is within a factor of the original Bayes Factor under 2n factor introduced due to the change of measure c can the natural probability measure. Hence, the algorithm 2n outweigh the gain made by the concentration of the divides the Bayes Factor by the factor h if the Bayes -nb probability measure e . This is not surprising because Factor is larger than one. If the Bayes Factor is less our construction thus far does not force the test not to than one, the algorithm multiplies the Bayes Factor by 2n bias against a sample in an intelligent way. That is, a the factor h . maliciously designed testing procedure could simply avoid the error prone regions of the design. To address Results and discussion this, we define the notion of a fair testing strategy that We applied our algorithm to two SDE models of tumor does not engage in such malicious sampling. dynamics from the literature. The first model is a single Definition 10. A testing strategy is h-fair (h ≥ 1) if dimensional stochastic differential equation for the and only if the geometric average of the implied Radon- influence of chemotherapy on cancer cells, and the sec- Nikodym derivatives over a number of samples is within ond model is a pair of SDEs that describe an immuno- a constant factor h of unity, i.e., genic tumor. Lefever and Garay model 1 p (X |u) μ i ≤ ≤ η Lefever and Garay [38] studied the growth of tumors η p (X |u) μ i i=1 under immune surveillance and chemotherapy using the following stochastic differential equation: Note that a fair test strategy does not need to sample from the underlying distribution in an i.i.d. manner. dx x βx x = r x (1 − ) − + x(1 − ) A cos (ωt) 0 0 However, it must guarantee that the probability of dt K 1+ x K observing the given behavior in a large number of + x (1 − ) W observations is not altered substantially by the non-i.i.d. sampling. Intuitively, we want to make sure that we bias Here, x is the amount of tumor cells, A cos(ωt) for each sample as many times as we bias against it. denotes the influence of a periodic chemotherapy treat- Our main result shows that such a long term neutrality ment, r is the linear per capita birth rate of cancer is sufficient to generate statistical guarantees on an cells, K is the carrying capacity of the environment, and otherwise non-i.i.d. testing procedure. b represents the influence of the immune cells. Note Definition 11. An h-fair test is said to be eventually 4 b that W is the standard Brownian Motion, and default fair if and only if 1 ≤ h <e ,where b is the constant in model parameters were those used in [38].. the exponential posterior concentration theorem. We demonstrate our algorithm on a simple property The notion of a eventually fair test corresponds to a of the model. Namely, starting with a tumor consisting testing strategy that is not malicious or adversarial, and of a billion cells, is there at least a 1% chance that the is making an honest attempt to sample from all the tumor could increase to one hundred billion cells under events in the long run. under immune surveillance and chemotherapy. The fol- lowing BLTL specification captures the behavioral speci- Algorithm fication: Finally, we present our Statistical Verification algo- rithm(SeeFigure2)interms of ageneric non-i.i.d. 10 11 Pr (F (x> 10 )) ≥0.01 testing procedure sampling with random “implied” change of measures. Our algorithm is relatively simple Figure 3 contrasts the number of samples needed to and generalizes our previous Bayesian Statistical verifi- decide whether the model satisfies the property using i.i. cation algorithm [8] to non-i.i.d. samples using change d. and non-i.i.d.. As expected, the number of samples of measures. The algorithm draws non-i.i.d. samples required increases linearly in the logarithm of the Bayes from the stochastic differential equation under ran- factor regardless of whether i.i.d.ornon-i.i.d. sampling domly chosen probability measures. The algorithm is used. However, non-i.i.d. sampling always requires ensures that the implied change of measure is bounded fewer samples than i.i.d. sampling. Moreover, the differ- so as to make the testing approach fair. The variable n ence between the number of samples increases with the denotes the number of samples obtained so far and x Bayes factor. That is, the lines are diverging. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 7 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 Figure 2 Non-i.i.d. Statistical Verification Algorithm. The figure illustrated the non-i.i.d. Bayesian model validation algorithm. The algorithm builds upon Girsanov’s theorem on change of measure and Bayesian model validation. Figure 3 Comparison of i.i.d. and non-i.i.d. sampling. Non-i.i.d. vs i.i.d. sampling based verification for the Lefever and Garay model. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 8 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 We note that there are circumstances when our algo- For this model, we considered the following property: rithm may require more samples than one based on i.i. starting from 0.1 units each of tumor and immune cells, d. sampling. This will happen when the property being is thereatleast a1%chancethatthe number of tumor tested has relatively high probability. For example, we cells could increase to 3.3 units. The property can be tested the property that the probability of eradicating encoded into the following BLTL specification: the tumor is at most 1%. The i.i.d. algorithm required 1374 samples to deny this possibility while the non-i.i.d. Pr (F (x > 3.3)) ≥0.01 algorithm required 1526 samples. Thus, our algorithm is best used to examine rare behaviors. Default model parameters were those used in [39]. Fig- ure 4 contrasts the number of samples needed to decide Nonlinear immunogenic tumor model whether the model satisfies the property using i.i.d.and The second model we analyze studies immunogenic non-i.i.d.. The same trends are observed as in the pre- tumor growth [39,40]. Unlike the previous model, the vious model. That is, our algorithm requires fewer sam- immunogenic tumor model explicitly tracks the ples than i.i.d. hypothesis testing, and that the difference dynamics of the immune cells (variable x) in response to between these methods increases with the Bayes factor. the tumor cells (variable y). The SDEs are as follows: We also considered the property that the number of tumor cells increases to 4.0 units. We evaluated whether dx(t)= (a − a x (t)+ a x (t)y(t)) dt 1 2 3 this property is true with probability at least 0.000005 +(b (x(t) − x )+ b (y(t) − y )) dW (t) 11 1 12 1 1 under a Bayes Factor of 100, 000. The i.i.d.sampling dy(t)= (b y(t)(1 − b y(t)) − x(t)y(t)) dt 1 2 algorithm did not produce an answer even after obser- +(b (x(t) − x )+ b (y(t) − y )) dW (t) 21 1 22 1 2 ving 10, 000 samples. The non-i.i.d. model validation algorithm answered affirmatively after observing 6, 736 The parameters x an y denote the stochastic equili- 1 1 samples. Once again, the real impact of the proposed brium point of the model. Briefly, the model assumes algorithm lies in uncovering rare behaviors and bound- that the amount of noise increases with the distance to ing their probability of occurrence. the equilibrium point. Figure 4 Comparison of i.i.d. and non-i.i.d. sampling. Non-i.i.d. vs i.i.d. Sampling based verification for the nonlinear Immunogenic tumor model. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 9 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 References Discussion 1. Faeder JR, Blinov ML, Goldstein B, Hlavacek WS: Rule-based modeling of Our results confirm that non-i.i.d. sampling reduces the biochemical networks. Complexity 2005, 10(4):22-41. number of samples required in the context of hypothesis 2. Haigh J: Stochastic Modelling for Systems Biology by D. J. Wilkinson. Journal Of The Royal Statistical Society Series A 2007, 170:261-261[http:// testing – when the property under consideration is rare. ideas.repec. org/a/bla/jorssa/v170y2007i1p261-261.html]. Moreover, the benefits of non-i.i.d. sampling increase 3. Iosifescu M, Tautu P, Iosifescu M: Stochastic processes and applications in with the rarity of the property, as confirmed by the biology and medicine [by] M. Iosifescu [and] P. Tautu Editura Academiei; Springer-Verlag, Bucuresti, New York; 1973. divergence of the lines in Figures 2 and 3. We note, 4. Twycross J, Band L, Bennett MJ, King J, Krasnogor N: Stochastic and however, that if the property isn’trarethenanon-i.i.d. deterministic multiscale models for systems biology: an auxin-transport sampling strategy will actually require a larger number case study. BMC Systems Biology 2010, 4(34)[http://www. biomedcentral. com/17520509/4/34/abstract]. of samples than an i.i.d. strategy. Thus, our algorithm is 5. Kwiatkowska M, Norman G, Parker D: Stochastic Model Checking. In only appropriate for investigating rare behaviors. Formal Methods for the Design of Computer, Communication and Software Systems: Performance Evaluation (SFM’07). Volume 4486. Springer;Bernardo M, Hillston J 2007:220-270, LNCS (Tutorial Volume). Conclusions 6. Kwiatkowska M, Norman G, Parker D: Advances and Challenges of We have introduced the first algorithm for verifying Probabilistic Model Checking. Proc 48th Annual Allerton Conference on properties of stochastic differential equations using Communication, Control and Computing IEEE Press; 2010. 7. Langmead C: Generalized Queries and Bayesian Statistical Model sequential hypothesis testing. Our technique combines Checking in Dynamic Bayesian Networks: Application to Personalized Bayesian statistical model checking and non-i.i.d.sam- Medicine. Proc of the 8th International Conference on Computational Systems pling and provides guarantees in terms of termination, Bioinformatics (CSB) 2009, 201-212. 8. Jha SK, Clarke EM, Langmead CJ, Legay A, Platzer A, Zuliani P: A Bayesian and the number of samples needed to achieve those Approach to Model Checking Biological Systems. In CMSB Volume 5688 of bounds. The method is most suitable when the behavior Lecture Notes in Computer Science. Springer;Degano P, Gorrieri R of interest is the exception and not the norm. 2009:218-234. 9. Wald A: Sequential Analysis New York: John Wiley and Son; 1947. The present paper only considers SDEs with indepen- 10. Lehmann EL, Romano JP: Testing statistical hypotheses. 3 edition. Springer dent Brownian noise. We believe that these results can Texts in Statistics, New York: Springer; 2005. be extended to handle SDEs with certain kinds of corre- 11. Edgar GA: Radon-Nikodym theorem. Duke Mathematical Journal 1975, 42:447-450. lated noise. Another interesting direction for future 12. Younes HLS, Simmons RG: Probabilistic Verification of Discrete Event work is the extension of these method to stochastic par- Systems Using Acceptance Sampling. In CAV, Volume 2404 of Lecture Notes tial differential equations, which are used to model spa- in Computer Science. Springer;Brinksma E, Larsen KG 2002:223-235. 13. Younes HLS, Kwiatkowska MZ, Norman G, Parker D: Numerical vs. tially inhomogeneous processes. Such analysis methods Statistical Probabilistic Model Checking: An Empirical Study. TACAS 2004, could be used, for example, to investigate properties 46-60. concerning spatial properties of tumors, the propagation 14. Lassaigne R, Peyronnet S: Approximate Verification of Probabilistic Systems. PAPM-PROBMIV 2002, 213-214. of electrical waves in cardiac tissue, or more generally, 15. Hérault T, Lassaigne R, Magniette F, Peyronnet S: Approximate Probabilistic to the diffusion processes observed in nature. Model Checking. Proc 5th International Conference on Verification, Model Checking and Abstract Interpretation (VMCAI’04), Volume 2937 of LNCS Springer; 2004. Acknowledgements 16. Sen K, Viswanathan M, Agha G: Statistical Model Checking of Black-Box The authors acknowledge the feedback received from the anonymous Probabilistic Systems. CAV 2004, 202-215. reviewers for the first IEEE Conference on Compuational Advances in Bio 17. Grosu R, Smolka S: Monte Carlo Model Checking. CAV 2005, 271-286. and Medical Sciences (ICCABS) 2011. 18. Gondi K, Patel Y, Sistla AP: Monitoring the Full Range of omega-Regular This article has been published as part of BMC Bioinformatics Volume 13 Properties of Stochastic Systems. Verification, Model Checking, and Abstract Supplement 5, 2012: Selected articles from the First IEEE International Interpretation, 10th International Conference, VMCAI 2009, Volume 5403 of Conference on Computational Advances in Bio and medical Sciences LNCS Springer; 2009, 105-119. (ICCABS 2011): Bioinformatics. The full contents of the supplement are 19. Jha SK, Langmead C, Ramesh S, Mohalik S: When to stop verification? available online at http://www.biomedcentral.com/bmcbioinformatics/ Statistical Trade-off between Costs and Expected Losses. Proceedings of supplements/13/S5. Design Automation and test In Europe (DATE) 2011, 1309-1314. 20. Clarke EM, Faeder JR, Langmead CJ, Harris LA, Jha SK, Legay A: Statistical Author details Model Checking in BioLab: Applications to the Automated Analysis of T- Electrical Engineering and Computer Science Department, University of Cell Receptor Signaling Pathway. CMSB 2008, 231-250. Central Florida, Orlando FL 32816 USA. Computer Science Department, 21. Langmead CJ, Jha SK: Predicting Protein Folding Kinetics Via Temporal Carnegie Mellon University, Pittsburgh PA 15213 USA. Lane Center for Logic Model Checking. In WABI, Volume 4645 of Lecture Notes in Computer Computational Biology, Carnegie Mellon University, Pittsburgh PA 15213 Science. Springer;Giancarlo R, Hannenhalli S 2007:252-264. USA. 22. Jha S, Langmead C: Synthesis and Infeasibility Analysis for Stochastic Models of Biochemical Systems using Statistical Model Checking and Authors’ contributions Abstraction Refinement.. SKJ and CJL contributed equally to all parts of the paper. 23. Jeffreys H: Theory of Probability Oxford University Press; 1939. 24. Gelman A, Carlin JB, Stern HS, Rubin DB: Bayesian Data Analysis London: Competing interests Chapman & Hall; 1995. The authors declare that they have no competing interests. 25. Jha S, Langmead C: Exploring Behaviors of SDE Models of Biological Systems using Change of Measures. Proc of the 1st IEEE International Published: 12 April 2012 Conference on Computational Advances in Bio and medical Sciences (ICCABS) 2011, 111-117. Jha and Langmead BMC Bioinformatics 2012, 13(Suppl 5):S8 Page 10 of 10 http://www.biomedcentral.com/1471-2105/13/S5/S8 26. Øksendal B: Stochastic Differential Equations: An Introduction with Applications (Universitext). 6 edition. Springer; 2003 [http://www.amazon. com/exec/obidos/redirect?tag=citeulike07-20\&path=ASIN/3540047581]. 27. Karatzas I, Shreve SE: Brownian Motion and Stochastic Calculus (Graduate Texts in Mathematics). 2 edition. Springer; 1991 [http://www.amazon.com/ exec/obidos/redirect?tag=citeulike07-20\&path=ASIN/0387976558]. 28. Wang MC, Uhlenbeck GE: On the Theory of the Brownian Motion II. Reviews of Modern Physics 1945, 17(2-3):323[http://dx.doi.org/10.1103/ RevModPhys. 17.323]. 29. Girsanov IV: On Transforming a Certain Class of Stochastic Processes by Absolutely Continuous Substitution of Measures. Theory of Probability and its Applications 1960, 5(3):285-301. 30. Pnueli A: The Temporal Logic of Programs. FOCS IEEE; 1977, 46-57. 31. Owicki SS, Lamport L: Proving Liveness Properties of Concurrent Programs. ACM Trans Program Lang Syst 1982, 4(3):455-495. 32. Finkbeiner B, Sipma H: Checking Finite Traces Using Alternating Automata. Formal Methods in System Design 2004, 24(2):101-127. 33. Harel D: Statecharts: A Visual Formalism for Complex Systems. Science of Computer Programming 1987, 8(3):231-274[http://citeseerx.ist.psu.edu/ viewdoc/summary?doi=10.1.1.20.4799]. 34. Ehrig H, Orejas F, Padberg J: Relevance, Integration and Classification of Specification Formalisms and Formal Specification Techniques.[http:// citeseerx.ist. psu.edu/viewdoc/summary?doi=10.1.1.42.7137]. 35. Berger J: Statistical Decision Theory and Bayesian Analysis Springer-Verlag; 36. Jeffreys H: Theory of probability/by Harold Jeffreys. 3 edition. Clarendon Press, Oxford; 1961. 37. Choi T, Ramamoorthi RV: Remarks on consistency of posterior distributions. ArXiv e-prints 2008. 38. Lefever R, RGaray S: Biomathematics and Cell Kinetics Elsevier, North-Hollan biomedical Press; 1978, chap. Local description of immune tumor rejection,:333. 39. Horhat R, Horhat R, Opris D: The simulation of a stochastic model for tumour-immune system. Proceedings of the 2nd WSEAS international conference on Biomedical electronics and biomedical informatics BEBI’09, Stevens Point, Wisconsin, USA: World Scientific and Engineering Academy and Society (WSEAS); 2009, 247-252[http://portal.acm.org/citation.cfm?id = 1946539.1946584]. 40. Kuznetsov V, Makalkin I, Taylor M, Perelson A: Nonlinear dynamics of immunogenic tumors: Parameter estimation and global bifurcation analysis. Bulletin of Mathematical Biology 1994, 56:295-321[http://dx. doi. org/10.1007/BF02460644], 10.1007/BF02460644. doi:10.1186/1471-2105-13-S5-S8 Cite this article as: Jha and Langmead: Exploring behaviors of stochastic differential equation models of biological systems using change of measures. BMC Bioinformatics 2012 13(Suppl 5):S8. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit

Journal

BMC Bioinformatics – Springer Journals

Published: Apr 12, 2012

There are no references for this article.