Maximizing Information Gain for the Characterization of Biomolecular Circuits
Tim Prangemeier, Christian Wildner, Maleen Hanst, Heinz Koeppl
MMaximizing Information Gain for the Characterization ofBiomolecular Circuits
Tim Prangemeier*, Christian Wildner*, Maleen Hanst, and Heinz Koeppl
Department of Electrical Engineering and Information Technology,Technische Universität Darmstadt, Germanytim.prangemeier;christian.wildner;maleen.hanst;[email protected]
ABSTRACT
Quantitatively predictive models of biomolecular circuits are im-portant tools for the design of synthetic biology and molecularcommunication circuits. The information content of typical time-lapse single-cell data for the inference of kinetic parameters is notonly limited by measurement uncertainty and intrinsic stochasticity,but also by the employed perturbations. Novel microfluidic devicesenable the synthesis of temporal chemical concentration profiles.The informativeness of a perturbation can be quantified based onmutual information. We propose an approximate method to per-form optimal experimental design of such perturbation profiles. Toestimate the mutual information we perform a multivariate log-normal approximation of the joint distribution over parameters andobservations and scan the design space using Metropolis-Hastingssampling. The method is demonstrated by finding optimal pertur-bation sequences for synthetic case studies on a gene expressionmodel with varying reporter characteristics.
KEYWORDS
Optimal experimental design, Synthetic biology, Molecular commu-nication, Chemical reaction networks, Information gain, Molecularprogramming, Information theory
ACM Reference Format:
Tim Prangemeier*, Christian Wildner*, Maleen Hanst, and Heinz Koeppl.2018. Maximizing Information Gain for the Characterization of Biomolec-ular Circuits. In
NANOCOM ’18: NANOCOM ’18: ACM The Fifth AnnualInternational Conference on Nanoscale Computing and Communication, Sep-tember 5–7, 2018, Reykjavik, Iceland.
ACM, New York, NY, USA, 6 pages.https://doi.org/10.1145/3233188.3233217
Synthetic biology enables the programming of living cells on themolecular level and promises to be the basis of a range of newtechnologies. Quantitatively predictive models of the involvedbiomolecular processes facilitate the design of synthetic parts andmolecular communications circuits. Assessment of the kinetic model * Both authors contributed equally.Permission to make digital or hard copies of all or part of this work for personal orclassroom use is granted without fee provided that copies are not made or distributedfor profit or commercial advantage and that copies bear this notice and the full citationon the first page. Copyrights for components of this work owned by others than theauthor(s) must be honored. Abstracting with credit is permitted. To copy otherwise, orrepublish, to post on servers or to redistribute to lists, requires prior specific permissionand/or a fee. Request permissions from [email protected].
NANOCOM ’18, September 5–7, 2018, Reykjavik, Iceland © 2018 Copyright held by the owner/author(s). Publication rights licensed to theAssociation for Computing Machinery.ACM ISBN 978-1-4503-5711-1/18/09...$15.00https://doi.org/10.1145/3233188.3233217 parameters is often characterized by a large degree of uncertainty,not only due to measurement uncertainty but also inherent cell-heterogeneity and intrinsic stochasticity of the involved process[5]. Continuous-time Markov chains provide a reasonable approxi-mation of the dynamics of stochastic reaction networks [16], suchas biomolecular circuits.The experimental process of perturbing and observing cells forthe purpose of parameter inference is shown schematically in Fig. 1.Time-lapse fluorescence readouts 𝑦 ( 𝑡 ) of single-cells contain infor-mation about the temporal behaviour of biomolecular circuits uponsome chemical perturbation 𝑢 ( 𝑡 ) [4, 21]. Bayesian inference canfor example be employed to calibrate model parameters [27]. Theinformation gained is not only limited by experimental constraints,measurement uncertainty and inherent heterogeneity within popu-lations of clonal cells, but also by the dynamics of the perturbation.In the past, technical constraints have limited the employed per-turbations to a single step or pulse. Recent advances in microfluidicsenable the generation of well defined temporal concentration pro-files [1, 22]. In this study, we exploit this new capability and findapproximately optimal perturbation profiles 𝑢 ( 𝑡 ) for the inferenceof biomolecular model parameters. flow u ( t ) Microfluidics Input DesignInferenceSingle Cell TracesMicroscopy
Perturbation: concentration profile
In-vivo system: molecular network (kinetic parameters to be inferred) in yeast cells
Observation: time-lapse fluorescent reporter abundance for each cell u ( t ) ty ( t ) y ( t ) Figure 1: Schematic of the time-lapse single-cell ex-periments upon which parameter inference for thebiomolecular circuit models considered here is based.
Bayesian optimal experiment design for the inference of stochas-tic reaction networks was studied theoretically in [18] for the caseof complete observations. These considerations were extended toincomplete and noisy observations in [25], where informativenesswas quantified using the expected log-determinant of the posteriorcovariance matrix. A Monte Carlo gradient descent method wasused to approximately solve the resulting optimization problem. a r X i v : . [ q - b i o . M N ] J a n ANOCOM ’18, September 5–7, 2018, Reykjavik, Iceland T. Prangemeier*, C. Wildner*, M. Hanst and H. Koeppl
Other approaches focus on the Fisher information matrix as ameasure of informativeness. In [9] the linear noise approximationis used to obtain a differential equation for the Fisher informa-tion matrix. The method proposed in [19, 20] combines a Gaussianapproximation of the Fisher information with moment-based infer-ence [26].Here, we follow the Bayesian approach and use the mutual in-formation as a measure of informativeness. In the context of cellbiology, such an approach is discussed in [13]. Based on an ordi-nary differential equation model of the chemical kinetics, mutualinformation was estimated using Monte Carlo simulations.
The observed behaviour of fluorescent reporter proteins stemsfrom the interaction of various biomolecules in a system of chem-ical reactions. We consider the dynamics of such a system as acontinuous-time Markov chain (CTMC) X ( 𝑡 ) on the time interval [ ,𝑇 ] , assuming the system is well mixed.The species 𝒳 , ..., 𝒳 𝑛 interact in 𝜈 reactions. Each reaction isassociated with a stoichiometric change vector v 𝑖 encoding thenet change of the molecule abundance and a propensity function ℎ 𝑖 ( x ) = 𝑐 𝑖 𝑔 𝑖 ( x ) . Here, 𝑐 𝑖 is the kinetic rate constant and 𝑔 𝑖 ( x ) is apossibly nonlinear function of the state. In this study, we considermass-action kinetics where 𝑔 𝑖 ( x ) essentially corresponds to theproduct of the involved reactant abundances. To simplify notation,we also introduce the vector of the rate constants c = ( 𝑐 , . . . , 𝑐 𝑣 ) 𝑇 .Intuitively, ℎ 𝑖 ( x ) 𝑑𝑡 is the probability that reaction 𝑖 fires in theinfinitesimal time interval 𝑑𝑡 . Consequently, the combined propen-sity of any reaction occurring follows as ℎ ( x ) = (cid:205) 𝜈𝑖 = ℎ 𝑖 ( x ) andthe waiting time between consecutive reaction events is exponen-tially distributed. In particular, if 𝜏 𝑗 denotes the time of the 𝑗 -threaction event, the waiting time is given by (( 𝜏 𝑗 − 𝜏 𝑗 − ) | X ( 𝜏 𝑗 − ) = x ) ∼ Exp ( ℎ ( x ) . The probability that the 𝑖 -th reaction has caused the event at 𝜏 𝑗 isgiven by 𝑝 ( X ( 𝜏 𝑗 ) = x + v 𝑖 | 𝑋 ( 𝜏 𝑗 − ) = x , 𝜏 𝑗 − 𝜏 𝑗 − ) = ℎ 𝑖 ( x ) ℎ ( x ) . The above equations lie at the heart of the stochastic simmulationalgorithm that allows to generate statistically exact sample pathsof the process [2]. More information on stochastic modelling ofchemical reactions systems can be found in [24].In practice, we cannot measure full sample paths. Instead, onetypically obtains noisy readouts at the sample times 𝑡 , . . . , 𝑡 𝑁 . Themeasurement 𝑌 ( 𝑡 𝑛 ) = 𝑦 𝑛 is connected to the latent state X ( 𝑡 𝑛 ) bythe observation model 𝑝 ( 𝑦 𝑛 | x ( 𝑡 𝑛 )) . We combine all measurementsinto a vector y = ( 𝑦 , . . . , 𝑦 𝑁 ) 𝑇 .The central goal of model calibration is to infer the kinetic rateparameters c from the noisy measurements y . We adopt a Bayesianapproach and choose a prior distribution 𝑝 ( c ) over the model pa-rameters. This prior distribution reflects the uncertainty about themodel parameters before the data y is obtained. Assuming thatthe measurements are independent given the latent path, the joint density over all quantities becomes 𝑝 ( c , x , y ) = 𝑝 ( c ) 𝑝 ( x | c ) 𝑁 (cid:214) 𝑛 = 𝑝 ( 𝑦 𝑛 | x ( 𝑡 𝑛 )) . where 𝑝 ( c ) refers to prior distribution over the parameters and x denotes the full latent sample path.In Bayesian inference, the quantity of interest is the parameterposterior 𝑝 ( c | y ) = ∫ 𝑝 ( c , x , y ) 𝑑 x 𝑝 ( y ) . This distribution is intractable because it requires the evaluation ofinfinite dimensional integrals with respect to the latent paths X .Several sampling-based approaches to approximate the aboveposterior have been developed recently [7, 27]. Since these ap-proaches are computationally intensive, here we focus on an opti-mal design approach that circumvents the posterior altogether.As a specific example that will also serve as the basis for severalcase studies, consider the following gene expression modelG off 𝑐 𝑢 ( 𝑡 ) −−−−− ⇀↽ −−−−− 𝑐 G on mRNA 𝑐 −−−−− ⇀ mRNA + P G on 𝑐 −−−−− ⇀ G on + mRNA P 𝑐 −−−−− ⇀ P mRNA 𝑐 −−−−− ⇀ ∅ P 𝑐 −−−−− ⇀ ∅ . Here, transcription is regulated by a single gene G that can switchbetween on and off states. As is common in synthetic biology, thetransition from the inactive to the active gene state is modulatedby an external perturbation 𝑢 ( 𝑡 ) . The input 𝑢 ( 𝑡 ) corresponds to theconcentration of the inducer molecule in the cell.The mRNA is transcribed from the gene and serves as a templateto produce an inactive protein P . An additional conversion isnecessary to obtain the active reporter protein P . This can, forexample, reflect maturation or folding times and enables us toconsider the properties of various reporter dynamics. The outcome Y of an experiment depends on the unknown modelparameters C and the perturbation 𝑢 ( 𝑡 ) . We express this in termsof a probabilistic model 𝑝 𝑢 ( 𝒚 | c ) . As illustrated in Fig. 2, we con-sider the inference task as an abstract communication problem.In this picture, the probabilistic model 𝑝 𝑢 ( 𝒚 | c ) corresponds toa noisy channel with input C , where the unknown rate constantchosen by nature from the prior distribution 𝑝 ( c ) . The inferencestep corresponds to decoding in the communication picture.From an information theoretic point of view, the mutual informa-tion 𝐼 𝑢 ( C ; Y ) quantifies the average information about C containedin the noisy measurement Y . Hence, a principled objective functionfor parameter inference is given by 𝐽 ( 𝑢 ) = 𝐼 𝑢 ( C ; Y ) = ∫ 𝑝 𝑢 ( c , y ) log 𝑝 𝑢 ( c , y ) 𝑝 ( c ) 𝑝 𝑢 ( y ) 𝑑 c 𝑑 y (1)and the optimal design choice 𝑢 ∗ is given by 𝑢 ∗ = arg max 𝑢 𝐽 ( 𝑢 ) [13]. The quantity 𝐽 ∗ = max 𝑢 𝐼 ( C ; Y 𝑢 ) resembles the capacity ofa noisy channel. However, to compute the channel capacity, onemaximizes with respect to the input distribution for a fixed channelmodel, whereas in optimal design, the roles are reversed. In the aximizing Information Gain for the Characterization of Biomolecular CircuitsNANOCOM ’18, September 5–7, 2018, Reykjavik, Iceland communication picture, we may see 𝐽 ∗ as the inference capacity ofthe experiment 𝑝 𝑢 ( y | c ) and 𝑢 ∗ as the capacity achieving design.In the classical approach to Bayesian experimental design ([14])the objective function 𝐽 ( 𝑢 ) is given as 𝐽 ( 𝑢 ) = ∫ 𝑝 𝑢 ( y ) (cid:20)∫ 𝑝 𝑢 ( c | y ) log 𝑝 𝑢 ( c | y ) 𝑝 ( c ) 𝑑 c (cid:21) 𝑑 y = E Y u [ 𝐷 𝐾𝐿 [ 𝑝 𝑢 ( c | y ) || 𝑝 ( c )]] , (2)where the Kullback-Leibler divergence between posterior and priordistribution corresponds to the reduction in uncertainty (or gain ininformation) by the particular observation Y = y . Since the outcome Y is unknown a priori, the expectation must be taken with respectto all possible outcomes. It is straightforward to show that Eqn. 2and 1 are equivalent. Experiment Y u Inference uC ˆC p u ( y | c ) Figure 2: Schematic depiction of the channel capacity approach tooptimal experimental design.
Approximations to optimal design are often based on Eqn. 2 andtarget the posterior 𝑝 ( c | y ) [23]. Here, we take an approach thatfocuses on Eqn. 1 and approximates the joint 𝑝 ( c , y ) . The idea behind this approximation is to choose a suitable para-metric family 𝑝 𝝓 ( c , y ) and find the closest representation of thejoint distribution by minimizing Kullback-Leibler divergence withrespect to the second argumentˆ 𝝓 = arg min 𝝓 𝐷 𝐾𝐿 [ 𝑝 ( c , y ) || 𝑝 𝝓 ( c , y )] . (3)This approach has favourable information theoretic properties inthe context of belief approximation [12].In the case of a stochastic chemical reaction system, a suitablechoice of the prior over the rate constants is provided by a productGamma distribution [27]. As an approximate model we choose amultivariate log-normal distribution. Let z = ( c , y ) denote the jointvector of parameters and observations with 𝑀 = 𝜈 + 𝑁 components.Then, the log-normal model corresponds to the density function 𝑝 𝝁 , 𝚺 ( z ) = | 𝜋 𝚺 | − 𝑀 (cid:214) 𝑘 = 𝑧 − 𝑘 exp (cid:20) − ( log ( z ) − 𝝁 ) 𝑇 𝚺 − ( log ( z ) − 𝝁 ) (cid:21) where the logarithms act component-wise on the vector z . Here, 𝝁 ∈ R 𝑀 and 𝚺 ∈ R 𝑀 × 𝑀 is a positive semidefinite matrix.The multivariate log-normal distribution has several desirableproperties in the present context. First, it is restricted to positivevalues. Second it has heavier tails than the multivariate normaldistribution. In contrast to a product Gamma distribution it alsoallows to model correlations between variables. Setting 𝝓 = ( 𝝁 , 𝚺 ) , a short calculation shows that the optimalparameters in the sense of Eqn. 3 are given byˆ 𝝁 = E [ log z ] , ˆ 𝚺 = E (cid:104) ( log z − ˆ 𝝁 ) ( log z − ˆ 𝝁 ) 𝑇 (cid:105) . Though the expectation with respect to the target distribution 𝑝 ( c , y ) cannot be evaluated analytically, it is straightforward toestimate the results using samples from the distribution 𝑝 ( c , y ) .Since the perturbation 𝑢 ( 𝑡 ) modulates one of the rate constants,a stochastic simulation algorithm capable of dealing with time-dependent propensities has to be used [2]. Alternatively, one canintroduce a dummy species with high firing rate ensuring regularupdates of the propensity in segments where 𝑢 ( 𝑡 ) is small. Wefollow the latter approach here.The log-normal distribution permits a closed form solution ofthe mutual information [11]. In particular, if 𝑝 ( c , y ) is multivariatelog-normal with parameters 𝚺 = (cid:18) 𝚺 𝑪𝑪 𝚺 𝑪𝒀 𝚺 𝒀 𝑪 𝚺 𝒀 𝒀 (cid:19) , the mutual information is given by 𝐼 ( C ; Y ) =
12 log (cid:18) | 𝚺 𝑪𝑪 || 𝚺 𝒀 𝒀 || 𝚺 | (cid:19) . (4)Since mutual information is invariant under one-to-one transfor-mations [10], we would obtain the same result by first applying alog transformation on C and Y and then performing a Gaussianapproximation in the log domain.The mutual information above may be rewritten as 𝐼 ( C ; Y ) =
12 log (cid:18) | 𝚺 𝑪𝑪 ||( 𝚺 C | Y | (cid:19) , where Σ C | Y = 𝚺 CC − 𝚺 CY 𝚺 − YY 𝚺 YC . This implies that changes in theposterior variance are exponentially related to the mutual informa-tion.Note, that since Σ CC is independent of 𝑢 , the approximate ob-jective function 𝐽 ( 𝑢 ) is equivalent to the classical D-optimalitycriterion in the log domain [6]. In our setup, the input 𝑢 ( 𝑡 ) modulates the activation rate of thegene. Experimentally, the perturbation 𝑢 ( 𝑡 ) is created by varyingthe concentration of a chemical signal over time. In practice, thepossible perturbations may be subject to various constraints due totechnical limitations or toxicity to the cells.In this study, perturbations are limited to piece-wise constantfunctions with fixed interval length and constant ℓ -norm. The timeinterval of the experiment [ ,𝑇 ] is divided into 𝐾 segments. Wedenote with 𝑢 𝑘 the input level at the 𝑘 -th sub-interval (Fig. 3). Thedesign space is parametrized as 𝒰 = (cid:40) ( 𝑢 , . . . , 𝑢 𝐾 ) 𝑇 ∈ R 𝐾 + : 𝐾 ∑︁ 𝑘 = 𝑢 𝑘 = 𝐸 (cid:41) where 𝐸 > u . ANOCOM ’18, September 5–7, 2018, Reykjavik, Iceland T. Prangemeier*, C. Wildner*, M. Hanst and H. Koeppl i npu t l eve l "step" segment 1 segment 2 segment K time in minutes i npu t l eve l "pulse" Figure 3: Schematic representation of typical ’step’ and ’pulse’ per-turbations, the triangles indicate reporter measurement timepoints;the input sequence consists of 𝐾 segments of given duration (in thiscase 20 min). Note, both depicted perturbations have | | u | | = . The restricted design space 𝒰 still contains uncountably manyperturbations. To explore 𝒰 , we propose a sampling algorithmbased on the Metropolis-Hastings method.We realize this by introducing an artificial probability density 𝑝 ( u ) ∝ exp ( 𝛽𝐼 u ( C , Y )) (5)on the design space where 𝛽 > u ( 𝑡 ) , we generate a trial value u ′ from a suitable proposal distribution 𝑞 ( u ′ | u ( 𝑡 ) ) . The proposedmove is accepted with probability 𝛼 given by 𝛼 = min (cid:40) , 𝑝 ( u ′ ) 𝑝 ( u ( 𝑡 ) ) 𝑞 ( u ( 𝑡 ) | u ′ ) 𝑞 ( u ′ | u ( 𝑡 ) ) (cid:41) . (6)For a fixed number of input segments 𝐾 , the design space 𝒰 cor-responds to the set of positive piece-wise constant functions withfixed step size and constant norm. This is equivalent to the standard 𝐾 -dimensional probability simplex. A suitable proposal distributionon this space is given by a Dirichlet distribution ([8]) of the form 𝑞 ( u ′ | u ( 𝑡 ) ) = DIR ( u ′ | 𝑏 u ( 𝑡 ) ) where 𝑏 is a tuning constant controlling the variance of the proposal.With this proposal and 𝑝 ( u ) as in Eqn. 5, the acceptance ratio 𝛼 canbe evaluated analytically. The resulting formula is straightforwardbut lengthy and thus not presented here. Metropolis-Hastings iteration I ( C , Y ) simulationmoving average Figure 4: Example run of the Metropolis-Hastings sampler with sim-ulated annealing showing the convergence to an approximate opti-mum of the mutual information. The thick line indicates a movingaverage with window size 500.
In an optimal design scenario, we are primarily interested inthe modes of 𝑝 ( u ) because by construction the maximum of 𝑝 ( u ) coincides with the maximum of 𝐼 u . Hence, simulated annealing ([15]) is employed by recursively increasing 𝛽 (Eqn. 5) during thesimulation to achieve the result shown in Fig. 4. The annealing shiftsprobability mass from flat regions near the tails towards the modeof the distribution until convergence to an approximate optimum). We apply the proposed method to the gene expression model de-scribed in Section 1.1 with the parameter configuration in Table 1(gene copy number of one). We consider not only the proposed bestperturbation sequences, but also the worst, to gain some intuitive in-sight into how these come to be. First, the reference case results arepresented, followed by the same model with a ’faster’ reporter. Theeffect of the ℓ -norm on the approximately optimal input vectors u ∗ is analysed. Finally, we analyse the effect of changing the reporterkinetics for the characteristic best found perturbation sequences. Table 1: Reference parameter configuration for the transcriptionmodel and coefficient of variation (CV) for priors.
Param. 𝑐 𝑐 𝑐 𝑐 𝑐 𝑐 𝑐 Mean / s − .
006 0 .
005 1 0 .
02 0 .
01 0 .
01 0 . − . .
005 0 . . − − We assume prior knowledge of the rate constants for all simu-lations (for example from a previous inference experiment). TheGamma distributed prior is given by 𝑝 ( c ) = 𝜈 (cid:214) 𝑖 = Γ ( 𝑎 𝑖 , 𝑏 𝑖 ) with 𝑎 𝑖 = CV − 𝑖 and 𝑏 𝑖 = 𝑐 𝑖 · CV 𝑖 ; 𝑐 𝑖 and the corresponding coeffi-cient of variation CV 𝑖 from the parameter configuration given inTable 1 (unless otherwise stated).The input is temporally discretized into 𝐾 = ℓ -norm is heldconstant for each optimization scenario. Simulated measurementsof the reporter abundance corrupted by independent log-normalobservation noise are recorded at intervals of 10 min, in line withcommon experimental protocols (for example [21]). Input sequencesare optimized with regard to the inference of parameters 𝑐 to 𝑐 .The joint log-normal model imposes strong assumptions on thedistribution 𝑝 ( c , y ) because it implicitly enforces a linear relation-ship in the log domain between the parameters and the measure-ment for a fixed perturbation u .Due to the complexity of the model, it is challenging to assessthe quality of such an approximation rigorously. As a simple sanitycheck, we considered histograms of marginal measurement distribu-tion 𝑝 ( 𝑦 𝑖 ) at individual time-steps and scatter plots of the pairwisemarginals 𝑝 ( 𝑐 𝑗 , 𝑦 𝑖 ) of the joint distribution. These are deemed tofollow the log-normal distribution well.A more rigorous validation of the proposed method, for exampleby quantitative comparison with the earlier approaches, is left forfuture work, however, the characteristics of the optimal perturba-tions found with our method are in agreement with those reportedin [25]. aximizing Information Gain for the Characterization of Biomolecular CircuitsNANOCOM ’18, September 5–7, 2018, Reykjavik, Iceland We first consider optimizing the input for the reference model,with a reporter that matures somewhat more quickly than thefastest known fluorescent proteins (1 / 𝑐 ≤ 𝑇 𝑚 , expected maturationtime 𝑇 𝑚 ≈ nu m b e r o f r e po r t e r m o l ec u l es time in min i npu t u * (t) optimal inputstep inputleast optimal input time in min optimal inputleast optimal input Figure 5: Mean reporter abundance (top) for the optimized per-turbation sequence u ∗ (bottom); reference (left) and (right) modelwith larger reporter maturation and decay rates ( 𝑐 = . − , 𝑐 = .
004 s − ). The best sequence is an early pulse, while a late pulse is foundto coerce the least informative behaviour of the cells. The mean re-porter trajectory in the optimal case exhibits dynamics that includeboth an accumulation and a decay phase, whereas the step perturba-tion and the least optimal variant merely capture reporter accumu-lation and omit the decay phase. These findings are in agreementwith results reported in [25]. The characteristics of the optimal andleast optimal perturbation sequences found are similar for slowerreporters with maturation times of up to 45 min, corresponding totypical in-vivo reporters in yeast.’Faster’ reporters are available, which follow a biomolecuar sys-tems dynamics without maturation delays. Translocation reportersshuttle mature fluorescent proteins between the nucleus and cy-toplasm of a cell within 1 min [3]. We simulate their behaviour byincreasing 𝑐 and 𝑐 associated with reporter maturation and decay.The characteristic optimal and least optimal inputs for a fasterreporter are shown in Fig. 5 (right) with the corresponding meanreporter trajectories. The faster reporter enables two accumulation-decay cycles to be observed within the allotted measurement du-ration. The optimal perturbation is a double pulse, with an earlyand a later pulse. The two reporter accumulation-decay cycles arepositioned as far apart as possible such that the reporter nearsextinction at the final measurement timepoint.In contrast to the result for the reference model, the least optimalperturbation for the faster reporter is not a late pulse but similar tothe standard step perturbation with slight fluctuations. The proteintraces for both the step and the late pulse perturbations achievea steady state by the second measurement timepoint and remainthere throughout the experiment. The effect of varying the ℓ -norm of the input sequence on the char-acter of optimal and least optimal perturbations was investigated.The results are depicted in Fig. 6. No substantial effect was foundfor the standard reporter, beyond a widening of the early pulse (Fig.6 A to C). At || u || / 𝐾 = Figure 6: Optimal input sequences and reporter trajectories withvarious ℓ -norms ( | | u | | / 𝐾 = , and , A to C respectively) for thereference model, and for the faster reporter model | | u | | / 𝐾 = (D). The character of the optimal perturbation is also conserved forstronger stimulation of the faster reporter (between Fig. 5 (right)and Fig. 6 D). An early pulse is followed by a decay phase and a late,somewhat wider, pulse, which absorbs most of the additional inputstrength. The least optimal perturbation, however, is markedlydifferent from the step found for || u || / 𝐾 = ℓ -norm for both models. In both cases the poor input sequencesachieve higher peaks than the optimal ones. For high peaks instimulation the gene is rarely turned off and therefore strong mod-ulation of the gene-on rate 𝑐 becomes irrelevant as the propensityis zero in the gene-on state. The optimization of the input perturbations in the previous sec-tions found two characteristic sequences: the ’early pulse’ and the’double pulse’. To gain an intuitive understanding of when to usewhich sequence and how experimental constraints imposed by thereporters influence the inference capacity, we now analyse the ef-fect of reporter kinetics on the mutual information 𝐼 ( C , Y ) . Theratio of reporter maturation and decay rates ( 𝑐 / 𝑐 = const . ) isheld constant. The reporters considered in the previous sections ANOCOM ’18, September 5–7, 2018, Reykjavik, Iceland T. Prangemeier*, C. Wildner*, M. Hanst and H. Koeppl have ratio of 𝑐 / 𝑐 = 𝑐 / 𝑐 =
10 respectively. The mutual
Figure 7: Mutual information 𝐼 ( C , Y ) of characteristic perturbationsover the ratio of reporter production and maturation rates ( 𝑐 / 𝑐 )for a constant ratio of maturation and decay rates; large ratio denotenominally ’faster’ reporters and all perturbations satisfy | | u | | = . information for the two characteristic perturbation profiles overvarious ratios of reporter production and maturation rates ( 𝑐 / 𝑐 )are shown in Fig. 7. The results indicate that for slower reportersthe early pulse is most informative. For faster reporters, however,the double pulse is the most informative.At the extreme ends of the reporter speed scale, both very slowand very fast reporters were found to deliver lower informationgains, seemingly regardless of perturbation strategy. The very slowreporters are hindered by the time constraint 𝑇 =
120 min, whilefaster reporters tend to decay toward extinction quickly.These results indicate that it may generally be beneficial to in-clude as many reporter production-decay cycles into an experimentas possible. However, optimization of the input perturbations isexpected to lead to the most informative experimental results. Theproposed method for maximizing the information gain for can bothdetermine whether to employ a single or a double pulse, or optimizethe input perturbation for a given reporter.
We propose an approximate method to perform optimal experi-mental design of chemical perturbation profiles for the inferenceof biomolecular circuit parameters fields of synthetic biology andmolecular communication. Our method scans the entire input pa-rameter space, finding both the most and least optimal input se-quences. This may be beneficial to gain deeper insights into whatconstitutes a good perturbation for the continuous-time Markovchains considered, as well as information on how intermediateprofiles perform and how wide the optimal peaks are.An advantage of our method is that it does not require an infer-ence step. Forward simulation of the process suffices to determinethe approximately optimal perturbation designs. Furthermore, theMarkov chain Monte Carlo framework employed can easily beextended to other scenarios and optimize other experimental pa-rameters such as the timing of measurements, for example.While we are cautious with regard to interpreting the results dueto the complexity of the system investigated, the results attainedfor two test cases showed two distinct characteristic perturbations.These were found to depend on the properties of the employedreporters (for a given system). Relatively fast reporters benefitedfrom repeated stimulation, while slower reporters achieved the most informative results for a strong early perturbation. In thefuture, our tool may be employed to generalize the characteristicsof ’optimal’ perturbations as the basis for intuitive design rules.
ACKNOWLEDGMENTS
This work is supported by LOEWE CompuGene.
REFERENCES [1] A. Ainla, I. Gözen, O. Orwar, and A. Jesorka. 2009. A microfluidic diluter basedon pulse width flow modulation.
Anal. Chem.
81, 13 (2009), 5549–5556.[2] D. F. Anderson. 2007. A modified next reaction method for simulating chemicalsystems with time dependent propensities and delays.
The Journal of ChemicalPhysics
Nat. Commun.
PLoS One
9, 6 (Jun 2014), e100042.[5] M. B. Elowitz, A. J. Levine, E. D. Siggia, and P. S. Swain. 2002. Stochastic geneexpression in a single cell.
Science (80-. ).
Model-oriented design of experiments . LectureNotes in Statistics, Vol. 125. Springer New York, New York, NY.[7] A. Golightly and D. J. Wilkinson. 2011. Bayesian parameter inference for sto-chastic biochemical network models using particle Markov chain Monte Carlo.
Interface Focus
1, 6 (2011), 807–820.[8] N. L. Johnson, S. Kotz, and N. Balakrishnan. 2000. Continuous multivariatedistributions : 1. Models and applications. (2000).[9] M. Komorowski, M. J. Costa, D. A. Rand, and M. P. H. Stumpf. 2011. Sensitivity,robustness, and identifiability in stochastic chemical kinetics models.
Proc. Natl.Acad. Sci. U. S. A.
Phys. Rev. E
69 (Jun 2004), 066138. Issue 6.[11] Kvalseth, T. 1982. Some informational properties of the lognormal distribution(Corresp.).
IEEE Transactions on Information Theory
28, 6 (Nov. 1982), 963–966.[12] R. H. Leike and T. A. Enßlin. 2017. Optimal belief approximation.
Entropy
19, 8(2017).[13] J. Liepe, S. Filippi, M. Komorowski, and M. P. H. Stumpf. 2013. Maximizing theinformation content of experiments in systems biology.
PLOS ComputationalBiology
9, 1 (Jan 2013), 1–13.[14] D. V. Lindley. 1956. On a measure of the information provided by an experiment.
Ann. Math. Stat.
27, 4 (1956), 986–1005.[15] J. S. Liu. 2001.
Monte Carlo strategies in scientific computing . New York.[16] D. A. McQuarrie. 1967. Stochastic approach to chemical kinetics.
J. Appl. Probab.
4, 03 (1967), 413–478.[17] R. Milo and R. Phillips. 2015.
Cell biology by the numbers . Taylor & Francis Ltd.[18] P. Nandy, M. Unger, C. Zechner, and H. Koeppl. 2012. Optimal perturbationsfor the identification of stochastic reaction dynamics.
IFAC Proc. Vol.
45 (2012),686–691.[19] J. Ruess and J. Lygeros. 2015. Moment-based methods for parameter inferenceand experiment design for stochastic biochemical reaction networks.
ACM Trans.Model. Comput. Simul.
25, 2, Article 8 (Feb. 2015), 25 pages.[20] J. Ruess, A. Milias-Argeitis, and J. Lygeros. 2013. Designing experiments tounderstand the variability in biochemical reaction networks.
Journal of The RoyalSociety Interface
10, 88 (2013).[21] C. Schneider, L. Bronstein, J. Diemer, H. Koeppl, and B. Suess. 2017. ROC’n’Ribo:Characterizing a Riboswitching Expression System by Modeling Single-Cell Data.
ACS Synth. Biol.
6, 7 (Jul 2017), 1211–1224.[22] M Unger, S S Lee, M Peter, and H Koeppl. 2011. Pulse width modulation ofliquid flows: towards dynamic control of cell microenvironments. (2011), 1567–9.[23] I. Verdinelli and K. Chaloner. 1995. Bayesian experimental design : a review.
Stat.Sci.
10, 3 (1995), 273–304.[24] D. J. Wilkinson. 2012.
Stochastic modelling for systems biology . CRC Press/Taylor& Francis. 335 pages.[25] C. Zechner, P. Nandy, M. Unger, and H. Koeppl. 2012. Optimal variationalperturbations for the inference of stochastic reaction dynamics. In . IEEE, 5336–5341.[26] C. Zechner, J. Ruess, P. Krenn, S. Pelet, M. Peter, J. Lygeros, and H. Koeppl.2012. Moment-based inference predicts bimodality in transient gene expression.
Proceedings of the National Academy of Sciences