Identifying Bayesian optimal experiments for uncertain biochemical pathway models

Pharmacodynamic (PD) models are mathematical models of cellular reaction networks that include drug mechanisms of action. These models are useful for studying predictive therapeutic outcomes of novel drug therapies in silico. However, PD models are known to possess significant uncertainty with respect to constituent parameter data, leading to uncertainty in the model predictions. Furthermore, experimental data to calibrate these models is often limited or unavailable for novel pathways. In this study, we present a Bayesian optimal experimental design approach for improving PD model prediction accuracy. We then apply our method using simulated experimental data to account for uncertainty in hypothetical laboratory measurements. This leads to a probabilistic prediction of drug performance and a quantitative measure of which prospective laboratory experiment will optimally reduce prediction uncertainty in the PD model. The methods proposed here provide a way forward for uncertainty quantification and guided experimental design for models of novel biological pathways.


Introduction
The intricate and lengthy process of discovering novel drugs stands to gain significant advantages through the application of computational methods. 1 These methods provide efficient and cost-effective avenues for exploring expansive chemical spaces, predicting molecular properties, and optimizing potential drug candidates.Dynamic models of biological systems (e.g., pharmacokinetic (PK) or pharmacodynamic (PD) models) are useful tools for studying biomolecular processes.Such models have been applied to study inter-and intra-cellular phenomena, including regulatory, metabolic, and signalling processes within human cells. 234These models are derived from first-principles approximations of the complex dynamics that occur in vivo.The utility of these biological models lies in their simplicity; tractably capturing qualitative system behaviors at the expense of prediction accuracy.Although these dynamic models are often used in preclinical drug development to determine in vivo drug response, 5 they are rarely applied in drug discovery.Biological pharmacodynamic models are valuable tools in both the optimal design and validation phases of identifying novel drug candidates, making them crucial for computational drug discovery.Furthermore, regulatory agencies are beginning to accept in silico studies as part of the validation and testing of new medical therapies and technologies. 67This underscores a need to formalise approaches for uncertainty quantification and model improvement for pharmacodynamic models.
Uncertainty quantification is a critical analysis that is often overlooked in computational drug design, despite its relevance and available software tools. 8Previous studies have applied methods of uncertainty quantification to models for calibrating algorithmic parameters in automated diabetes treatment systems, 9 a PK/PD cancer model for the antivascular endothelial growth factor in the present of a cancer therapeutic agent, 10 the calibration of dose-response relationships for the γ-H2AX assay to predict radiation exposure, 11 and peptide-MHC affinity prediction using a data-driven model. 12In the aforementioned studies, experimental or clinical data was readily available for the proposed uncertainty quantification procedures.For many PD models, the scarcity of quantitative data for calibrating these models has led to significant challenges in uncertainty quantification.This limited availability of experimental data results in non-unique and/or unconstrained parameter estimations, leading to issues of nonidentifiability. 13Researchers have acknowledged the presence of significant parameter correlations, nonidentifiability, and parameter sloppiness in dynamic biology models, where parameters can vary over orders of magnitude without significantly affecting model output. 14Given these limitations, it is evident that incorporating uncertainty quantification to biological dynamic models, even in the absence of available data, is crucial if they are to be used for computational drug discovery and evaluation.
Optimal experimental design (OED) has been proposed as a method for improving model parameter estimates in biological pharmacodynamic models.In OED, existing measured data is used to determine what new experiment(s) should be done to best reduce parameter uncertainties. 15Several works have applied different flavors of optimal experimental design to dynamic biological models.Work done by Liepe et al. 16 demonstrates OED to maximize expected information gained in a signaling pathway model featuring 6 differential equations, 4 uncertain parameters, and selecting between 5 experimental outcomes.Researchers in Bandara et al. 17 showed that a sequential experimental design for a cell signaling model, coupled with fluorescence microscopy-based experimental data, was able to reduce uncertainty in some parameters, while it increased in others due to nonidentifiability.This work also provides an example of selecting optimal experiments based on an experimentally relevant metric.This approach is in opposition to the more common purely information-theoretic approach.We note that the model used in this study is relatively small, only featuring 4 differential equations, 4 uncertain rate constants, and selecting between 3 experimental possibilities.Researchers in Eriksson et al. 18 utilized approximate Bayesian computation to identify approximate posterior distributions for uncertain model parameters.This work presents a larger scale problem, with 34 reactions for 25 species, 34 uncertain parameters, and 6 experimental possibilities from which to select.This approximate approach is applicable to scenarios where the the data-generation distribution (likelihood function) is not available in analytic form.
In this paper, we utilize methods from Bayesian optimal experimental design using exact Bayesian inference and Hamiltonian Monte Carlo (HMC) to recommend experiments for mechanistic model improvement.Our focus is on a dynamic model of programmed cell death, i.e., apoptosis.This model predicts synthetic lethality in cancer in the presence of a PARP1 inhibitor.Therefore, when we refer to selecting an optimal experimental design, we mean "which species should we measure from a PARP1-inhibited cell experiment to improve confidence in simulated lethality predictions of a given inhibitor?"Our goal is to first understand the impact of parameter uncertainty on the predictive reliability of this model.Then, the aim is to recommend experimental measurements that can be done in the lab to maximize confidence in simulated lethality predictions.To do this, we conduct parameter inference using HMC for large volumes of data generated from simulated experiments.This way, we can determine which experimental data, in expectation, would most reduce uncertainty in modeled drug performance predictions.Such an approach provides researchers with a quantitative method for linking physical experiment design to model improvement in the absence of existing data.
It is important to note that while our optimal experimental design study is conducted with simulated experimental data, and the results are thus subject to our modeling assumptions, the dynamic model underpinning our analysis contains biologically relevant and measurable parameters, thus offering insights that can guide future experimental endeavors.Furthermore, the framework we lay out here is general and applicable to any relevant biological modeling setting where experimental data is obtainable.
Formally, the contributions of this work are threefold: (1) Applying Bayesian optimal experimental design and highperformance computing (HPC) to a large-scale coupled ODE system of 23 equations and with 11 uncertain parameters and choosing between 5 experimental designs -the largest PD system considered thus far in a fully Bayesian OED framework of which the authors are aware; (2) Proposing novel decision-relevant metrics for quantifying the uncertainty in therapeutic performance of novel inhibitor drugs at a given concentration (IC 50 ); and (3) Identifying optimal experiments that minimize the uncertainty in therapeutic performance (i.e., simulated lethality) as a function of the inhibitor dosage of interest.
As we will show, there is preference for measuring activated caspases at low IC 50 to reduce uncertainty in the probability of cell death.At larger IC 50 , our results show that this uncertainty is maximally reduced by collecting data the mRNA-Bax concentrations.Therefore, we conclude that the decision as to which species to experimentally measure must be selected by considering trade-offs between the predicted BOED objectives and the inhibitor viability in vitro.The approach and results presented here thus bridge model-guided optimal design with uncertainty quantification, advancing drug discovery through uncertainty-aware decision-making.

Methods
The problem addressed in this work is to understand the impact of model parameter uncertainty on the predictive capacity of the PD model for PARP1-inhibited cellular apoptosis; a representative model for cancer drug evaluation.And given the effects of uncertainty and lack of experimental data, we wish reestablish a predictive model by recommending measurements that can be made in the lab to better constrain the model parameters.
The workflow to achieve the stated goals is as follows: 1. Acquire many (synthetic) experimental measurements for each prospective measurable species in the model 2. Construct prior probability distributions (i.e., beliefs regarding model parameter distributions in the absence of data) 3. Conduct parameter estimation via Bayesian inference using the model, data, and prior probabilities for each prospective experiment and over multiple data samples 4. Compute expected drug performance predicted by the model given posterior probability distributions (i.e., updated parameter distribution beliefs after incorporating data) 5. Rank and recommend experiments based on a metric that quantifies reliability in model predictions.
A graphical depiction comparing traditional parameter calibration to the proposed approach using simulated experimental data is shown in Figure 1.

Prior prediction Posterior prediction Posterior parameter distribution
Species "A" Prior vs. posterior inhibitor performance  Traditional model parameter calibration via Bayesian inference is shown in Figure 1a.The model-derived prior predictive distribution of concentration profiles for Species "A" is constrained by experimental data (in black).This leads to an updated probability density function for the model parameters θ referred to as the posterior distribution.The posterior distribution is then used to quantify inhibitor performance via probability of achieving apoptosis in simulation across an array of IC 50 values.We can then compare the expected inhibitor drug performance between the prior and posterior uncertainty.More accurate predictions by the posterior are thus expected because real data was used to update the model.

Concentration
The methodology utilized in the present work is presented in Figure 1b.Here, we do not have access to real, measured data to constrain our model.However, we can use the PD model and an expert-derived error model to generate an ensemble of simulated experimental data.In simulating and calibrating to a large volume of simulated experimental data, we account for uncertainty in what could be measured in reality (in black).This then leads to an ensemble of posterior distribution predictions for the parameters θ.We obtain ensembles of posterior parameter distributions for several potential measurable species in the model (e.g., Species "B" and Species "C").From these, we can compute an ensemble of inhibitor performances as measured by probability of apoptosis across a range of IC 50 values.This implies a distribution on expected inhibitor performances as predicted by "measuring" different species and using the simulated measurements for calibration.By comparing the distributions in expected inhibitor performance across different measured species, we can identify which experiment is "optimal" by selecting the one that leads to the greatest reduction in uncertainty.

Mathematical Model
The general form of a pharmakodynamic model is shown in Eqn. 1, which represents a set of coupled ordinary differential equations.
In the rate-based modeling framework used here, variables y(t) ∈ R n are state variables representing the time-varying concentrations of the n chemical species of interest within the cell signalling pathway.The parameters κ ∈ R m are the parameter data, which includes the kinetic rate constants for the reaction network.Initial concentrations of species are those at t 0 , specified by the value of the vector y 0 ∈ R n .The system parameters and time-varying species concentrations are related by the functions f : R n → R, indexed by ∀i ∈ {1, . . ., n}.Also, we separate the input parameters in Equation 1 into two vectors: u, the control variables, and κ the uncertain model parameters.The pre-specified input parameters u are any experimental controls set a priori.In the model considered here, this is the half maximal inhibitory concentration (IC 50 ) for a novel PARP1 inhibitor molecule.IC 50 is a commonly used metric in experimental biology for quantifying a molecule's ability to inhibit a specific biological function.

Bayesian Optimal Experimental Design
Bayesian optimal experimental design is a statistical framework that aims to identify the most informative experimental conditions or measurements to improve parameter estimation and model predictions in the presence of uncertainty. 19,20  involves iteratively selecting experiments 21,22 that maximize the expected information gain or other objectives, such as mean objective cost of uncertainty, 23,24 given prior knowledge, model uncertainty, and the available resources.The goal is to strategically design experiments that provide the most valuable information via uncertainty reduction in the calibrated model predictions.
In our application, we want to identify an experimental design ξ (i.e., a specific species we can measure from the cell in the presence of a PARP1 inhibitor) from the set of feasible designs, Ξ, which optimizes information gain regarding uncertain parameters, θ, given the outcome of that experiment, y.Let us denote the set of all uncertain parameters within the ODE system, including initial conditions and rate constants, as the set Θ.For tractability reasons, we may choose to only consider a subset of the parameters which are deemed most significant via global sensitivity analysis and refer to the vector of key, significant parameters as θ ∈ Θ.We also denote the vector of experimental designs being considered as ξ ∈ X ⊂ Ξ, and |X | = N ξ is the number of design decisions to choose from in optimization.
Here, the individual experiments represent a species to measure in an experiment representing the process modeled in the PARP1-inhibited cell apoptosis PD model.We assume prior distributions on the uncertain parameters p(θ), and a likelihood of the data conditional on the parameters and controls variables, p(y|θ, ξ, u).For a general objective function of interest, J (ξ; u), where ξ are the decision variables and u inputs to Equation 1, the Bayesian optimal experimental design task is described by Equation 2.
Traditionally, this objective J is the expected information gain (EIG), a quantitative measure of entropy reduction between the prior and posterior probability distributions of θ. 25,26 Under such a formulation, the objective function becomes a nested expectation computation under the marginal distribution of experimental outcomes.The posterior distributions are often then estimated using a Markov Chain Monte Carlo (MCMC) algorithm to compute the information gained.For this application, we are not interested in optimal experiments which simply minimize relative entropy in the prior and posterior distributions in θ, as is the standard practice in BOED.Instead, we wish to select optimal experimental designs that minimize uncertainty in lethality predictions made by the ODE model under different inhibitor drug IC 50 , i.e., different u ∈ U.In this model setting, lethality (i.e., probability of cell apoptosis) can be viewed as a probabilistic analog to the traditional dose-response curve and is defined in Equation 3.
To that end, we devise two metrics for quantifying uncertainty in cell lethality predictions at different IC 50 values.
The first is the uncertainty in the probability of triggering cell apoptosis at a given IC 50 , represented as σ apop (u).A graphical representation of σ apop (u) across IC 50 is shown in Figure 2a.This quantity defines a characteristic height of the distribution of the lethality curve predictions at a specific IC 50 of interest.Smaller σ apop values thus imply that there is higher confidence in the PD model predictions at the IC 50 they are computed at.One may choose to only consider the σ apop at very small IC 50 , since we are interested in drug molecules which are most effective at low doses.
The second metric provides a quantification of the uncertainty in the inhibitor drug dosage (i.e., IC 50 ) that achieves a given probability of triggering cell apoptosis, and is represented as σ IC50 (T).A graphical representation of of σ IC50 (T) across IC 50 and for T = 0.75 is shown in Figure 2b.This metric defines a characteristic width of the lethality distribution at a target probability of triggering cell death, T ∈ [0, 1].A smaller σ IC50 value at a given T implies higher confidence in which IC 50 value achieves a given probability of killing the cell.More detail on these metrics is discussed in Appendix 6.3.In addition, we discuss the multi-objective treatment of these two metrics, and the sensitivity of the optimal recommended experiment to that treatment, in Section 4.

Application
We consider a dynamic PARP1-inhibited cell apoptosis models for both wild-type human cells published in. 27We note that this pathway can be found in some tumorous cells, as well.Thus, for the present study, the process that this model represents is the inhibition of PARP1 as a therapy for cancer.By applying an inhibitor via a specified IC 50 , inhibition of PARP1 halts the DNA strand break repair process, reducing damage repair pathways and potentially leading to cell apoptosis.Further information on the model and the complete system of differential equations is available in the Equation S12.
The model consists of a set of 23 ordinary differential equations (ODE) wherein a subset of the 27 rate constant parameters are considered uncertain.This model takes as inputs (1) an IC 50 value and (2) initial conditions and (3) kinetic parameter values.In the present study, we only consider output at a specific time point, (i.e., the final time point) when performing inference.This means that the likelihood we supply for the data follows a Normal distribution in the model prediction at t = τ , where τ is the time limit for the numerical integrator used to solve the ODE system at a given IC 50 , and with a known σ that is computed with a 10% error relative to the model prediction.Information for the uncertain parameters considered in this study was first elicited from expert understanding of the system and relevant sources.Nominal parameter values and confidence interval information is shown in Table S1.Details regarding the construction of all prior distributions on uncertain parameters considered in the present study are provided in Table S2.In addition, because experimental data for species concentrations within a single cell undergoing PARP1-inhibited apoptosis are not readily available, we utilize simulated experimental data for our study that is generated via the procedure described in Appendix Section 6.1.

Objective function specification
As noted previously, the quantity-of-interest in assessing the performance of PARP1 inhibitor is the predicted lethality.The predicted lethalities are then used to compute our uncertainty metrics or objective functions, σ apop (u) and σ IC50 (T).
Here, we specify the method for computing lethality, i.e., probability of triggering cell death.It has been shown that a activated caspase molecule count in the cell above a threshold of 1,000 is sufficient to trigger cell apoptosis. 16Therefore, we can compute whether or not cell apoptosis has been irreversibly triggered by the introduction of a PARP1 inhibitor by computing y max k where k =Casp-act.Thus, lethality, or P (apoptosis) can be computed by sampling the prior (posterior) probability distributions of the uncertain parameters, θ i ∼ p(θ) ∀i ∈ [1, . . ., N s ], specifying a value for IC 50 , and integrating the system of ODEs.Given the ensemble of N s maximum-caspase concentrations, P (apoptosis) can be calculated as shown in Equation 3.
Here, N s is the Monte Carlo sample count.When computed at all IC 50 of interest, we retrieve a vector of P (apoptosis|u) predictions ∀u ∈ U, producing a lethality curve for a given posterior parameter estimate.

Incorporating biological constraints
In addition to the expert-elicited prior probability information in Table S1, there is other biologically and therapeutically relevant information that we can incorporate into Bayesian parameter estimation to further constrain our posterior predictions.First, we know that for PARP1 inhibitor drugs with sufficiently high IC 50 values (i.e., low efficacy), we expect the repair signal from PARP1 to win out over the programmed cell death induced by p53.This means that at high IC 50 , the DNA double strand break will be repaired and the cell will live.We also know the opposite to be true: for low IC 50 , programmed cell death will be induced and the cell will die.In mathematical terms, we can say:  4.
This prior knowledge regarding the effect of a drug on our quantity of interest does not directly translate to knowledge regarding the uncertain model parameters.However, because the quantity of interest is computed directly from an output of the model (i.e., concentration of activated caspases) we can still include this prior knowledge in Bayesian inference in three equivalent ways: in the likelihood, the prior, or through post facto filtering of the estimated posterior chain.For the present application, we choose the posterior filtering approach due to numerical challenges in the HMC sampling algorithm used to estimate the posterior distribution.
The posterior filtering approach used involves sampling a larger chain (in this study, 50,000 samples) for each simulated experimental data point.If a given parameter sample in the posterior chain leads to the violation of the constraints in Equation 4 when specified in the forward simulation model under the nominal IC 50 (i.e., Talazoparib, see Appendix Section 6.1), that sample is removed from the chain.

Results and Discussion
In this section, we provide the key findings of our study applying Bayesian optimal experimental design to the PARP1inhibited cell apoptosis model.First, we conduct a prior sensitivity analysis, follow by posterior estimation and computation of our optimal experimental design objectives.

Prior Sensitivity
We conducted a prior sensitivity analysis computed for a quantity of interest, P (apoptosis).The purpose of this study is to understand how much reduction in uncertainty is required to achieve desired lethality curve predictions.The result of this study is shown in Figure 3. Here, the P (apoptosis) is computed by (1) drawing 1,000 random samples from prior probability distributions for each uncertain parameter, then (2) solving the model in Equation S12 at those parameter sample values and (3) computing P (apoptosis) via Equation 3 for all IC 50 = [0.001,0.01, 0.032, 0.1, 0.32, 1.0, 3.2, 10.0, 100.0].To determine the model sensitivity against the priors, we scale σ for each parameter by a multiplicative factor to narrow the support of the prior distribution around its nominal mean.In Figure 3a, we neglect our prior constraints and just solve the model under the unconstrained prior samples.Under these conditions, we see that given the nominal prior uncertainty (i.e., 1.0σ), the effect of IC 50 on lethality vanishes (i.e., it is a more-or-less flat line).From the results, it appears that the model requires at minimum a reduction in prior uncertainty of two orders of magnitude (i.e., 0.01σ) before a predictive model is restored via sigmoidal curvature in the prior lethality curve.This implies that achieving a reduction in uncertainty spanning two orders of magnitude across all uncertain parameters is imperative to differentiate therapeutic performance between different inhibitors using the PD model.
In Figure 3b, we conduct the same sensitivity analysis while also applying the constraints in Equation 4to restrict the model-predicted lethality values using biological insights.Under these conditions, we obtain lethality curve predictions that have the desired sigmoidal shape, showing low therapeutic performance at large IC 50 and high therapeutic performance at small IC 50 .However, it is also evident that significant reduction in parameter uncertainty (0.0001σ) is required to achieve apoptotic switch behavior at low IC 50 . 16 is worth noting that this sensitivity analysis involved altering all prior σ values by the same amount, and the existence of parameter correlations could potentially yield alternative outcomes in terms of the necessary uncertainty reduction.However.these results underscores the need for substantial reduction in parameter uncertainty in PD models to re-establish predictive ability.

Posterior Estimation
We estimate posterior probability distributions for the uncertain parameters using the proposed Bayesian optimal experimental design approach.Details regarding the underlying algorithm are provided in Appendix Section 6.2.The summary statistics for (1) mean posterior parameter estimates, (2) effective sample size per parameter, and (3) filtered chain length are shown in Table S7-S16.The filtering procedure leads to non-uniform posterior chain lengths and lower effective sample sizes, as shown in the summary statistics.We note that only approximately 1% of samples in the posterior chains pass the constraints enforced via Equation 4. This is indicative of a known challenge in applying feasibility constraints after sampling: that the vast majority of the parameter space explored is ultimately infeasible, even for sufficiently large chains.Therefore, in future work, we aim to explore methods such as Bayesian constraint relaxation 28 to enforce relaxed constraints on our probabilistic model while maintaining numerical stability in standard Hamiltonian Monte Carlo sampling algorithms.Still, the average ESS per parameter is still close to the average chain length for each species, implying a high number of non-correlated samples in the posterior chains we computed.Additionally, we note that the assumption of Log-Normal prior probability distributions for all the model parameters in Table S1 is an assumption with ramifications on the final derived posterior distributions.

Value of Information Estimation
The ensemble of predicted lethality curves for each species in Table S3 are shown in Figure 4.It is clear from the plots that certain species induce less spread in lethality predictions across IC 50 (e.g., mRNA-Bax) while other exhibit more variability (e.g., Casp-pro).To quantitatively understand these trends, we compute the the values of σ apop (u) and σ IC50 (T).The plots in Figures 5a-5b show how these uncertainty metrics evolve over u and T, respectively.Tabulated results are available in Tables S4 and S5.It is seen in Figure 5a that there is a regime at small IC 50 (which is the regime with the greatest therapeutic viability) wherein the minimizer of σ apop changes.In particular, at IC 50 = 0.01 and IC 50 = 0.03, the ξ which minimizes the σ apop objective is Casp-act.At IC 50 > 0.03, the minimizer becomes mRNA-Bax for the rest of the IC 50 range considered.This can also be seen in Table S4, where the minimum value of σ apop in each column (i.e., for each IC 50 ) is demarcated with an * .
When considering the σ apop objective alone, and when only placing significance in reducing uncertainty at low IC 50 , the optimal decision ξ * = Casp-act.At IC 50 = 0.01, the potential reduction in uncertainty in P (apoptosis) under experimental designs (i.e., measurements) ξ * = Casp-act is 24% when compared to Bad-Bcl-xL, the species with the highest σ apop (0.01).We note again that reduction in σ apop directly translates to a reduction in the uncertainty regarding the predicted P (apoptosis) for a given inhibitor (i.e., at a given IC 50 ).This is analogous to reducing uncertainty in a given inhibitor drug's expected performance.In summary, this result implies that it is optimal to acquire and calibrate to experimental data regarding Casp-act concentrations if we wish to maximize confidence in the expected lethalities of inhibitors with with low IC 50 .
In the case of σ IC50 , as shown in Figure 5b, the trend of which ξ minimizes the objective remains consistent across threshold values T = [0.5, 0.6, 0.7, 0.8, 0.9].Thus, when considering the σ IC50 objective alone, and for any threshold, the optimal ξ * = mRNA-Bax.We note that reduction in σ IC50 directly translates to a reduction in the uncertainty regarding which IC 50 achieves the target P (apoptosis) threshold, T. This is analogous to reducing uncertainty in a which inhibitor drug achieves the desired expected performance.For a potent inhibitor with T = 0.9, selecting the optimal experimental design ξ * = mRNA-Bax over Casp-act (i.e., the species with the highest σ IC50 (0.9)) results in a substantial 57% reduction in uncertainty regarding the specific IC 50 associated with achieving this cell death probability.In summary, this result implies that it is optimal to acquire and calibrate to experimental data regarding mRNA-Bax concentrations if we wish to maximize model prediction confidence in which inhibitor dosage achieves a given expected lethality.
When combining these two metrics to determine an optimal experimental design (i.e., which protein to measure in the lab to reduce uncertainty in lethality curve predictions) we provide an optimal design ξ * against the following objectives: (1) σ IC50 at a threshold of 0.90, (2) σ apop at IC 50 = 0.01, and (3) weighted averages of both.In the case of the weighted averages, we take the objective to be defined as: w 1 σ IC50 + w 2 σ apop wherein w 1 + w 2 = 1.For the purposes of illustrating the sensitivity of the combined objective to the selection of weights w 1 , w 2 , we show this weighted average objective under w 1 = [0.01,0.1].These results are tabulated in Table S6.The results of the weighted average show that at w 1 = 0.01, Casp-act is optimal, while at w 1 = 0.1, mRNA-Bax is preferred.Although the weights used here are arbitrary, domain knowledge should be integrated in the assignment of weights, leading to optimal experiments that align with what is understood regarding the underlying biological pathways.
In each of the objectives considered, activated caspase or Bax mRNA are identified as the optimal experimental measurements that could maximally reduce model prediction uncertainty if used in calibration.Here, we provide some insight as to why this may be the case given the present study.A possible explanation for the switch from ξ * = Casp-act to ξ * = mRNA-Bax at low to intermediate IC 50 may be because of the time dynamics.If a higher IC 50 is needed to kill a cancer cell via apoptosis, then measuring a species that represent the precursor steps to triggering apoptosis may more informative than the activated caspase itself.It is also possible that mRNA-Bax performs best across a wider range of drug dosages because it is more proximal to the mechanism of action.Once PARP is inhibited, only a few reactions are needed to trigger mRNA-Bax formation, and ultimately, the activation of caspase.
These differing outcomes across the different experimental design objectives underscore the essential role of biological insights in differentiating among optimal experimental designs.

Conclusions
In this work, we utilized a realistic PD model for PARP1-inhibited cell apoptosis in a model-guided optimal experimental design workflow to reduce uncertainty in the performance of cancer drug candidates.We first showed that, under prior uncertainty, the ability for the model to discriminate between PARP1-inhibitor therapeutic performance disappears.Through extensive simulations and Bayesian inference, we have demonstrated the potential of our methodology to reduce parameter uncertainties and reestablish a discriminatory PD model for drug discovery.Additionally, our exploration of the PARP1 inhibited cell apoptosis model has yielded novel insights into species-specific dependencies with PARP1 inhibitor IC 50 , informing optimal measurement choices.Introducing therapeutic performance metrics based on the half maximal inhibitory concentration (IC 50 ), we provide a comprehensive evaluation of drug efficacy under parameter uncertainty.Specifically, we identified optimal experiments under three distinct objectives: (1) σ IC50 at a threshold of 0.90, (2) σ apop at IC50=0.01, and (3) a weighted average of both metrics.At low IC 50 , our results indicate that measuring activated caspases is optimal, with the potential to significantly reduce uncertainty in model-predicted cell death probability.And at high IC 50 , Bax mRNA emerges as the favored contender for experimental measurement.We note that our algorithmic approach, which includes collecting a large number of posterior samples and then rejecting parameter samples that violate biological constraints, leads to a majority of samples being rejected.Therefore, in future work, we will explore alternatives for enforcing biological constraints in Bayesian inference, such as relaxed constraints in the data likelihood.Additionally, we would like to identify lab experiments that can be conducted to acquire data which is compatible with our model so that we can validate the results of our OED and further constrain the model parameters.In summary, our combined metric-based approach not only elucidates the importance of tailored species measurements to minimize uncertainty in lethality curve predictions but also introduces a flexible framework for optimal experimental design using novel biological pathway models.Thus, our contributions bridge the critical gap between model-guided optimal design and uncertainty quantification, advancing the realm of computational drug discovery.

Simulated Experimental Data
Simulated experimental data was generated for each of the N ξ observable proteins by taking random parameter samples from their respective prior predictive distributions (see Figure S8) and modifying these samples via a relative error of 10% across all parameters.The value of 10% was selected as it is the typical standard deviations found in experimental measurements in the laboratory.First, N d random draws from the parameter prior distributions are made independently for each uncertain parameter, i.e.
Next, the θ sim i are used to solve the forward model in Equation S6 at an IC 50 =0.005083µM , which is the measured IC 50 for WT cells under the PARP1-inhibitor Talazoparib, as shown in Figure 6.This generates a "true" y.
Figure 6: The HCT116 dose response curves.Details on the study can be found in. 27e noisy y sim k,i ∀k = [1, . . ., N ξ ], ∀i = [1, . . ., N d ] used in inference in this work is generated by applying a multiplicative error ϵ, where we assume that the error is log-normally distributed.Therefore, Where the standard deviation σ is 0.1 for all species k.It follows that: To avoid numerical issues at 0 concentrations (i.e., y k,i = 0), we employ a parameter a defined such that a = b − min{y k,i }, where b is a small tolerance (in this study, b = 0.001).Leading to the final form of our simulated experimental data as:

Algorithmic approach
The workflow diagram depicted in Figure 7 outlines the systematic procedure employed for conducting Bayesian Optimal Experimental Design (BOED) within the context of the PARP1-inhibited cell apoptosis model.The process begins by specifying prior uncertainty descriptions for the key model parameters (p(θ) ∀θ ∈ Θ) determined by global sensitivity analysis results presented in. 27This subset is shown in Table S1.Next, synthetic experimental data (y sim k,i ) are generated for each experimental design (indexed by k = {1, . . ., N ξ }) under consideration, as listed in Table S3, and for i = {1, . . ., N d } samples from the prior distributions and a known error model (see Table S7).We note that we select the experimental designs listed in Table S3 because these species exhibited high prior predictive uncertainty, as seen in Figure S8.These synthesized data are used to constrain the uncertain model parameters via Bayesian Inference, where the posterior parameter probability distributions (p(θ|y sim k,i , ξ k )) are estimated using a HMC algorithm (i.e., a type of MCMC algorithm for sampling from an unknown distribution) for each k = {1, . . ., N ξ } and for each datum i = {1, . . ., N d }.Next, constraints based on biological feasibility are then applied to reject samples in the Markov chains that violate the expected behavior of the cell under extreme PARP1 inhibitor dosages.This leads to modified posterior estimates (p ′ (θ|y sim k,i , ξ k )).These posterior distributions are then utilized to compute the Quantity of Interest (QOI), specifically the probability of cell apoptosis (P (apoptosis|u)), which is assessed for various inhibitor IC 50 values (u ∈ U).These QOI values serve as the basis for computing the uncertainty metrics: σ apop (u) and σ IC 50k (T) for each experimental design.Finally, the uncertainty metrics are combined to make optimal experimental decisions, thereby providing recommendations for the most effective experimental designs.This comprehensive workflow guides the process of leveraging BOED to enhance our understanding of how we can optimally reduce uncertainty in the model predictions made by the PARP1-inhibited cell apoptosis model under different drug dosage conditions.

Uncertainty Metrics
The calculated lethality serves as the basis for deriving our fundamental experimental design objectives: the σ IC50 (T) and σ apop (u) uncertainty metrics.These metrics offer a quantifiable means of assessing the reduction in uncertainty associated with the various experimental designs, ξ.
The σ IC50 (T) metric is determined by fitting a 4-parameter sigmoidal curve, specifically the form P (apoptosis|IC 50 ) b i , to the computed P (apoptosis|u) i values for each IC 50 value and for each i ∈ [1, . . ., N d ].Subsequently, the intersection points of the sigmoidal fits with pre-specified P (apoptosis) thresholds (e.g., T = [0.5, 0.6, 0.7, 0.8, 0.9]) determine the associated IC 50 values.The σIC 50 (T) metric is then defined as the standard deviation over IC 50 values at which the sigmoidal fits intersect the threshold.
The σ apop metric is computed by taking the standard deviation over all predicted probabilities of apoptosis at a given IC 50 , i.e., P (apoptosis|u) i , ∀i = [1, . . ., N d ], ∀u ∈ U.This is defined in Equation 5, where P (apoptosis|u) is the arithmetic mean of the predicted probabilities of cell death at the specified inhibitor IC 50 , or u.

Implementation
The entirety of the workflow presented here was implemented using the programming language Julia. 32Specifically, we utilize the probabilistic programming library Turing.jl 29to implement our probabilistic model for Bayesian inference.
To sample the posterior chain and obtain posterior probability distributions, we use the No U-Turn Sampler (NUTS), a variant of HMC, with an acceptance rate of 0.65.The ODE model evaluations within the probabilistic model are conducted via the DifferentialEquations.jllibrary. 30e implementation of BOED we propose here exploits several elements of high-performance computing approaches to accelerate the optimization process.First, the implementation of the ODE model was translated into the low-level programming language, Julia, from the rule-based modeling language, BioNetGen. 31When performing ODE model simulations at the nominal parameter settings, the BioNetGen implementation requires 0.10 seconds of CPU time.In contrast, solving the ODE model under the same nominal parameter settings using Julia takes only 0.0018 seconds of CPU time.Consequently, our Julia-based model implementation exhibits an impressive 98% reduction in CPU time compared to the BioNetGen implementation.Although these individual execution times may already appear low, the acceleration of individual model evaluations is critically important in the context of Bayesian optimal experimental design.This is because these model evaluations are used in estimating the data likelihood across various parameter values, and as the number of these evaluations increases, the cumulative effect on computational time of Bayesian inference becomes more pronounced.Secondly, we accelerate the computation of each of the 100 posterior estimates per experimental design by using 100 CPU cores in parallel.This means that we can compute 100 posteriors simultaneously, decreasing the time-to-solution over a sequential approach.Thus, we show that the combined utilization of these high-performance computing strategies enhances the efficiency and practical applicability of our proposed Bayesian optimal experimental design framework within the context of a complex PD model.All instances of posterior sampling through NUTS were solved on the Institutional Cluster at the Scientific Data and Computing Center at Brookhaven National Library.Runs were done in parallel using 3 nodes at a time.

Experimental Designs
The list of measurable species considered as potential "experimental designs" in the present study are shown in Table 3 9 Tabulated Objective Data σ apop (0.01) σ IC50 (0.9)

Posterior Chain Summary Statistics
Each table in this section presents summary statistics for uncertain parameter values and effective sample size for each uncertain parameter inferred in Bayesian inference against a specific species measurement.Averages are computed over 100 synthetic experimental datum for that protein generated according to the procedure described in the STAR Methods.An additional table is provided per each species to list the min, max, and average chain size from all 100 posteriors used in obtaining these statistics.

Mathematical Model
This section contains the full system of differential equations that define the PARP1-inhibited cell apoptosis model as presented in. 27In the notation in Equation 6, ẋi refers to the time-varying rate of change of species i, where i = {1, . . ., 23}, x i is the species concentration in units of molec/cell, and p i , i = {1, . . ., 27} are the ODE parameters related to rates.See Tables 17 and 18 for a mapping from x i and p i notation to biological symbols and values.This mathematical model was generated from the rule-based modeling software, BioNetGen, using the previously published PARP1-inhibited cell apoptosis model in.18: Mapping between symbol and parameter name in the ODE model in Equation 6.There is no nominal value specified for p 23 =IC 50 because it is a decision variable and its value is subject to change within our studies.
Model parameter calibration with simulated experimental data.

Figure 1 :
Figure 1: Workflow for Bayesian parameter calibration using (a) real, measured data and (b) experimental data simulated from the model.

Figure 2 :
Figure 2: Graphical depictions of uncertainty metrics for probabilistic lethality curves in (a) σ apop (u) at u = 0.001 and (b) σ IC50 (T) at target probability threshold T = 0.75.The shaded green region represents a distribution of lethality curves that is computed from an ensemble of posterior distributions in the uncertain parameters θ.The mean of this distribution in lethality curves is shown in the square data points.

Figure 3 :
Figure 3: The probability of cell death P (apoptosis) plotted against IC 50 for different multiplicative factors of the standard deviation σ in the prior distributions (a) without and (b) with applying the constraints in Equation 4.

Figure 4 :
Figure 4: Ensemble of lethality curves predicted for each measurable species considered at different IC 50 values.Each individual curve represents a different lethality prediction obtained from a single synthetic experimental measurement.The mean over all 100 measurements is shown in black.

Figure 5 :
Figure5: Uncertainty metrics for synthetic lethality computed at different IC 50 using posterior distributions inferred using concentration data for the measurable species considered.

Figure 7 :
Figure 7: Workflow utilized for Bayesian experimental design for the PARP1-inhibited cell apoptosis model.

Figure 8 :
Figure 8: Prior predictive plots from 500 simulations sampled randomly from the prior distributions of the uncertain parameters and run through the forward ODE model to create concentration profile distributions.Units of concentration are in molecules-per-cell.

Table 1 :
Table of uncertain parameters and 90% confidence limits for the PARP1 inhibitor ODE model.All prior probability distributions are assumed to be log-normal to prohibit negative rate parameter values.Derived prior probability distributions for each uncertain parameter considered in the model are described in Table2.All priors are defined as Log-Normal probability distributions.For all parameters, the nominal value was taken as the mean.Because the lower and upper bounds on the uncertain parameters are not centered at the nominal value, the width of this CI range was used to compute the geometric distance from the nominal to new lower and upper bounds in log-space.

Table 2 :
Table of distributional information (mean, standard deviation) for Log-Normal priors on all uncertain parameters considered in this work.

Table 3 :
27Table of measurable species in the PARP1 inhibited ODE model.As in,27all initial values for these species is 0. The ODE model computes concentrations in dimensionless units [Molec.]representing molecules-per-cell.

Table 4 :
Table of σ apop results for each measureable protein in the PARP1-inhibited cell apoptosis model.The minimum value in each column is denoted with an * .

Table 5 :
Table of σ IC50 results for each measureable protein in the PARP1-inhibited cell apoptosis model.The minimum value in each column is denoted with an * .

Table 6 :
Table of optimal experimental designs, ξ * for each objective considered in the present study.

Table 7 :
Average expected parameter values and effective sample size from filtered chains for posteriors obtained via mRNA-Bax data.

Table 8 :
Chain length statistics (min, max, and average) after sample rejection for mRNA-Bax.

Table 9 :
Average expected parameter values and effective sample size from filtered chains for posteriors obtained via Bad-Bcl-xL data.

Table 10 :
Chain length statistics (min, max, and average) after sample rejection for Bad-Bcl-xL.

Table 11 :
Average expected parameter values and effective sample size from filtered chains for posteriors obtained via Casp-pro data.

Table 12 :
Chain length statistics (min, max, and average) after sample rejection for Casp-pro.

Table 13 :
Average expected parameter values and effective sample size from filtered chains for posteriors obtained via Bax-Bcl-xL data.

Table 14 :
Chain length statistics (min, max, and average) after sample rejection for Bax-Bcl-xL.

Table 15 :
Average expected parameter values and effective sample size from filtered chains for posteriors obtained via Casp-act data.

Table 16 :
Chain length statistics (min, max, and average) after sample rejection for Casp-act.