An in-silico analysis of experimental designs to study ventricular function: A focus on the right ventricle.

In-vivo studies of pulmonary vascular disease and pulmonary hypertension (PH) have provided key insight into the progression of right ventricular (RV) dysfunction. Additional in-silico experiments using multiscale computational models have provided further details into biventricular mechanics and hemodynamic function in the presence of PH, yet few have assessed whether model parameters are practically identifiable prior to data collection. Moreover, none have used modeling to devise synergistic experimental designs. To address this knowledge gap, we conduct a practical identifiability analysis of a multiscale cardiovascular model across four simulated experimental designs. We determine a set of parameters using a combination of Morris screening and local sensitivity analysis, and test for practical identifiability using profile likelihood-based confidence intervals. We employ Markov chain Monte Carlo (MCMC) techniques to quantify parameter and model forecast uncertainty in the presence of noise corrupted data. Our results show that model calibration to only RV pressure suffers from practical identifiability issues and suffers from large forecast uncertainty in output space. In contrast, parameter and model forecast uncertainty is substantially reduced once additional left ventricular (LV) pressure and volume data is included. A comparison between single point systolic and diastolic LV data and continuous, time-dependent LV pressure-volume data reveals that at least some quantitative data from both ventricles should be included for future experimental studies.

Introduction

Computational modeling, combined with invasive or non-invasive measurements, can forecast both the onset and worsening of cardiovascular disease [1-3]. More recently, multiscale models that account for cardiovascular physiology across multiple spatial scales have been developed [4,5]. The synergistic combination of in-vivo and in-silico methods have had notable success in understanding the progression of right ventricular (RV) failure in pulmonary hypertension (PH) [6-9]. The left ventricle (LV) and septal wall (S) are highly coupled to RV function [10]; hence, an impaired RV reduces biventricular energy efficiency and overall LV function [11]. The use of mechanistic models and their physiologically based parameters can reveal additional details of PH progression, especially when combined with highly informative in-vivo data. However, these computational models suffer from numerous parameters and limited, noisy data available for parameter inference and model calibration [12].

In these situations, a formal identifiability analysis can reveal which parameters to infer, and which data collection protocols are most informative for the model. There are two main types of identifiability. Parameters are considered structurally identifiable if the model output is unique for every unique parameter set. In addition to structural identifiability, parameters can also be practically identifiable if they can be uniquely determined from limited and/or noisy data. Structural identifiability assesses the model’s structure, and is determined using algebraic manipulations of the model [13-15] or by inferring parameters using noise-free, model generated data [1,16]. Parameters that are deemed structurally identifiable can be assessed for practical identifiability in the presence of noisy and limited data. This type of analysis is imperative to inform in-vivo experimental designs for the frequency or quality of measurements.

Several authors have considered parameter identifiability in the context of cardiovascular modeling [1,15,17,18]. Pironet et al. [15] pursued a structural identifiability analysis on a six-compartment model of the cardiovascular system. The study concluded that a combination of pressure and volume data was necessary to eliminate structural non-identifiability for the 13 parameters in their model. A follow up investigation by Pironet et al. used local sensitivity and profile likelihood analyses to conclude that only subset of parameters were practically identifiable from swine data in the vena cava, aorta, and LV [17]. The studies by Colunga et al. [18] and Harrod et al. [1] used models including the LV, RV, and both the systemic and the pulmonary circulations. The former [18] used local sensitivity analysis and Markov chain Monte Carlo (MCMC) methods to deduce identifiable parameter subsets given limited data from heart transplant patients. The latter study [1] utilized similar sensitivity and MCMC techniques, and tested for structural identifiability by examining the marginal posterior distributions for each parameter after fitting the model to noise-free, model generated data from patients with PH due to left heart failure. Both studies found practical identifiability issues in the full parameter set, and instead deduced a smaller subset of model parameters that were both identifiable and physiologically meaningful.

These prior studies did not consider a multiscale model with biventricular interaction. This latter component is especially important during the progression of PH and during chronic RV pressure overload [6,10]. The cutting edge reduced order model of biventricular interaction is the three-segment (“TriSeg”) model developed by Lumens et al. [8]. Two recent studies by van Osta [5,19] applied sensitivity analysis and uncertainty quantification methods to the TriSeg model, and identified which parameters were influential on model forecasts of RV, LV, and S wall strain. These investigations utilized non-invasive clinical data, whereas only a few studies have used the TriSeg model with in-vivo animal data [4,9,20,21]. Animal models of PH provide novel insight into PH progression [22,23], yet it is unclear how informative in-vivo data from these experiments are for calibrating computational models.

To address these gaps in knowledge, this study investigates practical parameter identifiability for a multiscale model of biventricular interaction and cardiovascular dynamics. We utilize sensitivity analyses, the profile likelihood, and MCMC techniques to deduce practical identifiability of the model. We focus on data obtained from four experimental designs; three that are common for monitoring animal models of PH and focus on the RV [23-25], and an additional design that utilizes dynamic pressure-volume data in both the LV and RV [26]. We generate both noise-free and noisy data from the model to test for practical parameter identifiability and analyze the output uncertainty in model simulations and several biomarkers of PH progression.

Materials and methods

Mathematical model

We consider a multiscale cardiovascular model describing sarcomere-level dynamics, biventricular interaction, and zero-dimensional (0D) hemodynamics. We summarize the mathematical model here and relegate individual component details to S1 Text.

The model consists of nine compartments: the systemic and pulmonary arteries and veins, the left and right atria, and a model accounting for interactions between the LV, RV, and S. A model schematic is provided in Fig 1.

Fig 1

Model schematic.

The computational model here consists of a lower order simulator of sarcomere dynamics within the left atrium (LA), left ventricle (LV), right atrium (RA), right ventricle (RV), and septum (S). The LV, RV, and S are simulated using the TriSeg model [8], and account for biventricular interaction. Lastly, a circuit model is used to describe the systemic arteries (SA) and veins (SV), as well as the pulmonary arteries (PA) and veins (PV).

Sarcomere model

The sarcomeres in the atrial, ventricular, and septal walls are modeled as two passive elastic elements in parallel with an elastic and contractile element in series [27]. The contractile sarcomere length, L (μm), and contractility, Γ (dimensionless), are dictated by ordinary differential equations [28]. As described by Lumens et al. [8], changes in sarcomere length, L (μm), are dependent on myocardial strain, ε (dimensionless), while changes in Γ depend on L and time t (s). Cardiac contractility is modeled as the sum of a rise and decay function, describing the binding of crossbridges, calcium fluctuations, and detachment of crossbridges during diastole. Active stress, G (KPa), is determined as a function of L and Γ, whereas passive stress due to structural properties of the extracellular matrix, G (KPa), and the giant protein Titin, G (KPa), are strictly a function of sarcomere length. The total stress generated from the sarcomere is then the sum of the active and passive stresses

This subcomponent of the model constitutes a total of 27 parameters: 13 shared between the two atria, 13 shared between the LV, RV, and S, and a parameter describing the time delay of atrial contraction (see S1 Text).

TriSeg Model

The sarcomere model is embedded within a cardiac tissue model of atrial dynamics and biventricular interaction (the “TriSeg” model [8]), and relates changes in blood volume V(t) (μl) to myocardial strain ε, using

Here, A (mm2) is the current mid-wall area of the chamber, A (mm2) is the reference mid-wall area, and z (dimensionless) is a curvature variable related to the ratio of wall volume, V (mm2), and radius of mid-wall curvature C (mm-1) [8]. Once ε has been calculated and the corresponding G is obtained from the sarcomere model, the mid-wall tension can be calculated as

A balance in axial and radial tensions, T and T (see S1 Text), is enforced providing two differential algebraic equations [29]. The cavity tensions are used to calculate the cavity pressures (see S1 Text). In total, the cardiac chambers and TriSeg model contribute two algebraic constraints in Eq (4), five wall volume parameters (V), and five reference area parameters (A).

Hemodynamics model

The systemic and pulmonary arteries and veins are modeled as compliant compartments, with resistance elements between each compartment or cardiac chamber [18,30]. In brief, changes in V, flow q, μl/s, and pressure p (KPa) are related via an electric circuit analogy where the subscripts in and out denote the compartments before and after a model component, V (μl) is the unstressed volume assumed at zero pressure (see S1 Text), C (μl KPa-1) is the vascular compliance, and R (KPa s μl-1) is the resistance between compartments. Finally, we model the two atrioventricular valves (mitral and tricuspid), the two semilunar valves (aortic and pulmonic), and the resistor between the systemic veins and right atrium as diodes

The hemodynamics model consists of eight differential equations for V(t), eight resistance parameters, and four compliance parameters.

Summary

The multiscale model consists of 18 differential equations (describing L, Γ, and V), two algebraic constraints (Eq (4)), and a total of 49 parameters. Due to the algebraic constraints, the model constitutes a system of differential algebraic equations (DAEs) and is solved using the variable-step, variable-order solver available in MATLAB (Mathworks; Natick, MA). Fig 2 shows nominal model predictions as well as the noise-corrupted data used in Bayesian parameter inference, discussed later.

Fig 2

Nominal simulations and noise corrupted data.

The nominal simulations are generated to match the data range reported by Philip et al. [23] in sham mice. Noise corrupted data is generated by adding additive, Gaussian errors with mean zero and a variance of 1.

Model sensitivity

Sensitivity analysis is an a posteriori identifiability method for determining which parameters are influential on a model output [15]. Local sensitivity analysis perturbs parameters one at a time, and typically utilizes finite-difference approximations [30,31]. In contrast, global sensitivity analysis samples parameters throughout the feasible parameter space, and includes variance based methods and screening methods [31,32]. We utilize a Morris screening analysis [33] in combination with a local, derivative based sensitivity analysis to determine parameter identifiability. We utilize Morris screening over variance based methods since the model parameter space is large (∈ℝ49). Prior studies have shown agreement between Morris’ indices and the total Sobol’ index [34], hence screening can be used to fix non-influential parameters.

The local sensitivity of a model output f with respect to a parameter, θ, is approximated by the centered difference where i = 1,2,…49 is the parameter index, f(t; ) is the quantity of interest from the model, Δθ is the step change in parameter value, and is the i-th unit vector. For time-dependent outputs, we consider the 2-norm of the model output, i.e. . We account for differences in parameter magnitude by computing the log-scaled parameter sensitivity [31,35]

The Morris’ screening approach computes the “elementary effects” where δ(ℓ) is the parameter step size describing the “levels” of effects. Choosing ℓ to be even provides a more symmetric sampling distribution [33], hence we choose ℓ = 60 giving δ≈0.51. Note that EE is a coarser approximation of model sensitivity than , but is qntified over a larger parameter space. We scale parameters from their original value to the interval [0,1] as done previously [34], and utilize the algorithm provided by Smith [12] to construct our sampling methodology. The indices from the Morris method are determined from K random initializations of the parameter vectors and are defined by

Here, μ is the average of EE, is an improved metric for average model sensitivity [34], and is the variance of EE. We use the combined index, , to measure a parameter’s influence [36].

Small values of either the local sensitivity index or the screening index M indicate that a parameter is non-influential, i.e. it has minimal effect on f. As discussed next, these indices assess whether a model parameter is practically identifiable.

Practical parameter identifiability

In this work, we assess practical parameter identifiability using three techniques. The first is through the local and global sensitivity metrics discussed above. Next, we consider the profile likelihood, which provides information about whether each θ is identifiable from a given set of data. Lastly, we use MCMC methods for Bayesian inference, and utilize the marginal posterior distributions to assess parameter identifiability.

Sensitivity based identifiability

Parameters that have little effect on the model output are considered practically non-identifiable, since they do not affect the quantity of interest [12], and should be fixed before conducting inference. We employ a two-part parameter fixing methodology using the results from Morris screening and local sensitivity analysis.

A parameter is deemed non-influential for all outputs f if its index M is less than the average for all parameters i = 1,2,…,P where f is one of the model outputs [5,32]. Parameters that are less than this threshold for all outputs are considered non-influential for inference and are fixed.

After using the Morris screening approach, the subset is analyzed by conducting a local sensitivity analysis around the nominal parameter values. The Fisher information matrix, , must be non-singular for gradient based parameter estimation, hence its utility in parameter identifiability [37]. If is invertible but has a large condition number (e.g., on the order of 1e8), then some of the sensitivities are nearly linearly dependent and the subset requires further reduction. We use an eigenvalue-eigenvector analysis method via the singular value decomposition (SVD) to determine which parameters cause the ill-conditioning of [14,38], and fix these parameters at their nominal value.

Profile likelihood

The most common and robust technique for assessing practical identifiability is the profile likelihood [13,15]. This technique increments a fixed parameter, θ, while minimizing the negative log-likelihood for all other parameters in the subset, i.e. Where y is the k-th data source, f is the corresponding model output, LL(y|) is the log-likelihood, is the noise variance for the data source, and N is the number of data points. The corresponding profile likelihood confidence intervals for θ are [13]

Each CI(θ) is constructed around the optimal estimate, *, and depends on the inverse cumulative distribution function of the chi-squared distribution, , with one-degree of freedom and confidence level α [13]. If PL(θ) is completely flat (e.g., CI(θ) is infinite), then θ is deemed structurally non-identifiable and cannot be uniquely determined due to model structure. If only one side of PL(θ) is flat, then θ is considered practically non-identifiable, and could become identifiable if more data was available for inference [16].

Bayesian inference

We assess the parameter identifiability in the presence of noise using Bayesian parameter inference. Using MCMC for Bayesian inference is more computationally expensive than gradient based optimization, but provides detailed insight into parameter relationships and avoids local minima in the likelihood [39-41]. We use the DRAM algorithm [42], which is described in depth elsewhere [31,43]. In short, the goal of MCMC is to approximate the posterior distribution where is the prior distribution, L(y|) is the likelihood, and the denominator of Eq (16) is a normalization factor. Model parameters are sampled from a proposal distribution to compute the likelihood L(y), where * is the proposed parameter values. The proposed parameter vector is accepted if the ratio of the likelihood values between * and the previous value of are greater than some random realization from a unit normal distribution. To reduce parameter stagnation or random-walk behavior, a second proposal parameter set is generated from a narrower distribution if * is rejected [43]. The DRAM algorithm updates the covariance matrix of the proposal after sequential adaption intervals, improving the proposed values of * [43].

We utilize DRAM on a set of noisy data, generated by the model at the nominal parameter values and corrupted with noise. To ensure adequate parameter space coverage and test the robustness of the MCMC, we first generate twelve random samples of our parameter subset and initialize a gradient based optimization that minimizes the residual sum of squared errors for the given experimental conditions (defined in the next section). Each optimal parameter vector, , is used as a starting value for an instance of DRAM, and the Hessian matrix obtained from the optimization is used as the initial covariance matrix to preserve possible sampling asymmetry [12]. We implement this using the freely available DRAM package developed Haario et al. [42] in MATLAB. In situations where the model is unstable or crashes, we return a large value for the residual sum of squares [39]. We assess parameter identifiability by visualizing the marginal posterior densities ); longer, unbounded tails in the posterior suggest issues with parameter identifiability. We assess MCMC convergence by looking at the median acceptance rate across chains as well as the potential scale reduction factor (PSRF) and multivariate PSRF (MPSRF). As suggested by Roy [44] we use a PSRF and MPSRF cutoff of 1.1 as an indicator of MCMC convergence.

Simulated experiments and additional outputs

Several experimental designs are commonly used for in-vivo PH studies [23,24]. We are interested in using the computational model to infer parameters indicative of heart function; hence, we consider different assortments of ventricular pressure and volume data. Our experimental designs are:

Dynamic measurements of RV pressure ();

Dynamic measurements of RV pressure and volume ();

RV pressure and volume measurements, as well as systolic and diastolic pressure and volume in the LV (); and

Dynamic measurements in both the RV and LV ().

The first two scenarios correspond to in-vivo recordings from pressure [24] or pressure-volume catheters [3,23]. The third includes additional information on the LV obtained by echocardiography [23]. Finally, the fourth experimental design represents a realistic, but underutilized, scenario that includes pressure-volume measurements in both the RV and LV [11,26].

We perform all sensitivity and identifiability analyses with respect to the pressure and volume forecasts considered in the four experimental designs above. Noisy pressure and volume data are generated by adding zero mean, white Gaussian noise, with a variance of 1 mmHg and 1 μl, respectively. Parameter subsets for each experimental design are contrasted, with a common subset determined across all designs. To better understand the consequences of limited data, we construct profile likelihood confidence intervals and analyze the parameter posterior distributions for each design. In the latter case, we compare the maximum a posteriori estimates with the known, data generating parameters. Lastly, we propagate uncertainties in the model parameters to simulated outputs via the posterior distributions. This includes LV, RV, and S engineering strain [19] as well as mean pulmonary artery pressure, RV stroke volume (the difference between end-diastolic and end-systolic volumes), pulmonary arterial elastance (the difference in mean pulmonary artery and left atrial pressure over stroke volume), RV end-systolic elastance, and RV ventricular-vascular coupling [23]. A graphical summary of the proposed parameter reduction workflow using sensitivity analyses, profile likelihood, and MCMC is provided in Fig 3. The source code of the mathematical model and relevant analyses can be found at https://github.com/mjcolebank/Colebank_Identifiability_2022.

Fig 3

Workflow schematic.

The initial set of 49 parameters is reduced to 38 due to apriori parameter fixing. Morris screening is used to confirm which parameters are on average the most influential on the four model outputs. This reduces the parameter set from 38 to 17 parameters. A local sensitivity analysis using the different experimental designs as the quantities of interest is used to determine if any parameters show local interdependence in their sensitivities, which suggests possible practical non-identifiability. If the Fisher information matrix constructed from the model sensitivity is ill-condition, the least influential parameters of the subset are fixed. This reduces the parameter subset from 17 parameters to a set of 13 parameters. Lastly, the parameters and experimental designs are subjected to profile likelihood analysis and MCMC to test for practical identifiability.

Results

Before beginning the sensitivity analysis, several parameters were excluded for physiological reasons. For instance, the reference length of the sarcomeres, L were excluded from analysis since these values are consistent across experimental designs. Table 1 summarizes the model parameters that are considered in our analyses. A detailed description of how parameter values are calculated can be found in the S1 Text. The ODE solver error tolerance is set to 10−12 to ensure smooth solutions, and the model is run for 60 cardiac cycles establish convergence to steady state. Simulation time ranges from 7–9 seconds depending on the parameters specified.

Table 1

Parameters, their description, and information regarding the sensitivity analyses.

Parameters with the subscript j have atrial and ventricular components.

Parameter	Description	Used in sensitivity analyses	Parameter	Description	Used in sensitivity analyses
V LA,wall	LA wall volume	X	τ offset,A	Offset of atrial systole
V LV,wall	LV wall volume	X	σ¯_act,j	Active stress scaling	X
V RA,wall	RA wall volume	X	σ¯_pas,j	Passive stress scaling	X
V RV,wall	RV wall volume	X	L s,pas,ref,j	Reference length for passive wall constituents
V S,wall	S wall volume	X	β pas,j	Stiffness of passive element	X
A m,ref,LA	LA reference area	X	k 1,j	Nonlinear scaling of Titin stiffness	X
A m,ref,LV	LV reference area	X	R a,val	Aortic valve resistance	X
A m,ref,RA	RA reference area	X	R m,val	Mitral valve resistance	X
A m,ref,RV	RV reference area	X	R p,val	Pulmonic valve resistance	X
A m,ref,S	S reference area	X	R t,val	Tricuspid valve resistance	X
L s,ref,j	Reference sarcomere length at zero strain		R vc	Vena Cava resistance	X
L s,iso,j	Elastic series element length in isometric state		R pv	Pulmonary venous resistance	X
v 0,j	Velocity of sarcomere shortening	X	R sys	Systemic circulation resistance	X
L sc,0,j	Contractile element length		R pulm	Pulmonary circulation resistance	X
Γrest,j	Resting contractility		C sa	Compliance of systemic arteries	X
τ rise,j	Rise in contractility scaling	X	C sv	Compliance of systemic veins	X
τ decay,j	Decay in contractility scaling	X	C pa	Compliance of pulmonary arteries	X
τ sys,j	Length of systole	X	C pv	Compliance of pulmonary veins	X

We ran the Morris screening algorithm using 100 randomized initializations. Fig 4 shows the parameter ranking M using the mean effect μ* and corresponding variance s2 for the RV and LV pressures and volumes. See S1 Text for individual results from the Morris screening analysis as well as parameter bounds for sampling. Sensitivity results were analyzed by comparing the parameter ranking M to the mean effect for each ventricular pressure and volume. All four compliances were consistently ranked within the most influential parameters, while other parameters describing cardiac chamber dynamics (e.g., A) varied with the output. We fixed parameters that were less influential than for all four outputs (i.e., pressure and volume in the RV and LV). This reduced our parameter subset from 38 to 17 parameters, shown in Table 2.

Fig 4

Sensitivity results from the Morris screening algorithm.

Parameter ranking is based on the index . The dotted line in each plot denotes the average model sensitivity for each output.

Table 2

Parameters deemed influential by Morris screening and included in the final subset after using a local sensitivity based practical identifiability analysis.

	Deemed important by Morris				Final subset
Parameter	p rv	p lv	V rv	V lv
V la,wall
V lv,wall		X		X	X
V ra,wall
V rv,wall			X		X
V s,wall			X		X
A m,ref,la	X				X
A m,ref,lv	X	X		X	X
A m,ref,ra
A m,ref,rv	X		X		X
A m,ref,s
v max,A
τ rise,A
τ decay,A
τ sys,A
σ¯_act,A
σ¯_pas,A
β pas,A
k 1,A
v max,V
τ rise,V		X			X
τ decay,V	X	X			X
τ sys,V	X	X	X	X	X
σ¯_act,V		X	X	X	X
σ¯_pas,V
β pas,V
k 1,V
R a,val
R m,val		X		X
R p,val
R t,val
R vc
R pv
R sys	X	X	X	X	X
R pulm	X		X		X
C sa	X	X	X	X
C sv	X	X	X	X
C pa	X	X	X	X	X
C pv	X	X	X	X

We conducted a local sensitivity analysis on the reduced subset of 17 parameters using the designs 1, 2, 3, and 4 as the quantity of interest. The local sensitivity of these designs with respect to the 17 parameters are used to construct the Fisher information matrix, . Using the SVD decomposition, we reduced the parameter subset until cond()≤108 for each design, providing a subset of 13 parameters deemed practically identifiable for all four designs. Parameters fixed by the SVD method included mitral valve resistance, R, and compliance in the systemic arteries, systemic veins, and pulmonary veins (C, C, and C, respectively). This final subset, shown in Table 2, was used in the profile likelihood and MCMC analysis.

Profile likelihood-based confidence intervals are constructed using the noise-free, model generated data. We construct the confidence intervals ±50% away from the true parameter value, with the confidence level cutoff for each design calculated using Eq (15) with an α = 0.95 confidence level. The profile likelihood results, displayed in Fig 5, show that only the last experimental designs, 4, provided finite confidence bounds for all 13 parameters. Sharp edges in the profile likelihood correspond to local minima and/or incompatible parameter sets corresponding to a failure in the DAE solver. For example, large values of V in combination with small values of A can lead to negative chamber volumes because of inadequate filling, making the parameter choice non-physiological. The parameters A, τ, τ, τ, R, and C were identifiable for all four experimental designs. The remaining seven parameters (V, V, V, A, A, σ, and R) varied in their identifiability with each experimental design.

Fig 5

Profile likelihood confidence intervals.

Confidence intervals are constructed by fixing one parameter and inferring all others over a range of values. Each row corresponds to a different experimental design. Note that the minimally informative experimental designs (1 and 2) have non-identifiable parameters, indicated by infinite or one-sided confidence bounds. In contrast, inclusion of LV data (3 and 4) remedy the issue of non-identifiable parameters in the set. Large deviations in the profile likelihood correspond to local minima and parameter sets that are incompatible for the system of DAE’s.

Noise corrupted data generated by the model is used in the likelihood defined in Eq (14). We use minimally informative priors (i.e., with a large variance) for each parameter and initialize the DRAM algorithm using the optimal parameter vector and estimated covariance matrix from twelve randomly selected initial guesses. MCMC is run for 50,000 iterations, with the initial 10,000 being left out as a “burn-in” period. We separate the results from MCMC into three groups: parameters representing the heart chambers’ geometry (Fig 6), parameters within the sarcomere model (Fig 7), and hemodynamic parameters in the circulatory model (Fig 8). Three of the twelve MCMC chains as well as the posterior distribution calculated using kernel density estimation are shown. The posterior distributions are relatively wide when only using RV pressure data (1), but additional data in the subsequent experimental designs reduce the posterior widths. All the marginal posterior distributions contain the true, data generating parameters, though some of the posteriors’ modes are unaligned with the true parameters. Additional pairwise plots, provided in S2 Text, suggest some correlation between variables. One chain using 3 shows tight, narrow correlations, likely due to a poor initialization during the optimization and inadequate exploration of the parameter space. However, the other eleven instances suggest minimal correlations between variables for 3. In general, 50,000 iterations appear sufficient for most of the MCMC results; however, the addition of static LV data with 3 causes some suboptimal mixing for the parameter σ. The PSRF for each parameter and the MPSRF for each design are provided in Table 3. These results suggest that 50,000 iterations of MCMC do not satisfy the cutoff of 1.1 as commonly used. Running an additional 25,000 iterations for each chain (results not shown) reduced the MPSRF slightly, but not below 1.1. The MCMC chains appear to converge quicker when using the most detailed experimental design, 4. The median acceptance rates for the twelve chains are 29.8%, 21.7%, 36.1%, and 49.0% for designs 1, 2, 3 and 4, respectively.

Fig 6

Chain iterations and marginal posteriors after MCMC for the TriSeg parameters.

The model parameters indicative of the TriSeg geometry (wall volume, V, and reference mid-wall area, A) are shown for each experimental design, corresponding to each column. The true, data generating parameters corresponding to the outputs in Fig 2 are shown as red lines. Three of the twelve initializations of MCMC are shown in different shades of gray. The marginal posterior distributions for the simplest experimental design (1) are much wider than the subsequent more informed experimental designs, suggesting an improvement in practical identifiability.

Fig 7

Chain iterations and marginal posteriors after MCMC for the sarcomere parameters.

Similar to Fig 6, three of the twelve MCMC instances are provided for the sarcomere parameters important for the rise, decay, and length of fiber shortening (τ, τ, and τ, respectively), and maximal active force generation (σ). Note that all four experimental designs (given by each column) provide sufficient information to the likelihood so that the true data generating parameters (in red) are within the marginal posteriors.

Fig 8

Chain iterations and marginal posteriors after MCMC for the hemodynamic compartment parameters.

As in Figs 6 and 7, three of the twelve MCMC instances are provided for systemic vascular resistance, R, pulmonary vascular resistance, R, and pulmonary arterial compliance, C. Though the marginal posteriors do contain the true parameters (in red) within the marginal posteriors for the simplest design (1, first column), additional data in the other designs substantially reduce posterior uncertainty.

Table 3

Potential scale reduction factor (PSRF) and multivariate PSRF (MPSRF) values calculated for each parameter and the 12 initializations of MCMC, respectively.

Parameter	PSRF f1	PSRF f2	PSRF f3	PSRF f4
V wall,LV	1.07	1.32	1.27	1.39
V wall,RV	1.28	1.22	1.41	1.22
V wall,S	1.13	1.08	1.07	1.04
A m,ref,LA	1.26	1.14	1.83	1.01
A m,ref,LV	1.16	1.14	1.08	1.19
A m,ref,RV	1.15	1.09	1.17	1.03
τ rise,V	1.10	1.05	1.05	1.22
τ decay,V	1.13	1.06	1.05	1.09
τ sys,V	1.07	1.11	1.13	1.23
σ act,V	1.12	1.18	1.87	1.35
R sys	1.10	1.17	2.71	1.05
R pulm	1.16	1.05	1.52	1.02
C pa	1.02	1.04	1.10	1.02
MPSRF	1.46	1.43	4.25	1.34

V

V

V

A

A

A

τ

τ

τ

σ

R

R

C

We propagate the uncertainties in model parameters to the outputs by subsampling from the posterior distributions. To account for any across chain variation, we draw fifty samples from the twelve different MCMC instances, giving 600 realizations from the posteriors. Fig 9 displays the noise-corrupted data, average response from the agglomerated samples, and one standard deviation from the average response. The results from the initial design, 1, show little uncertainty in RV pressure, but large uncertainty in forecasts of LV pressure and both chamber volumes. In contrast, 2, 3, and 4 show reduced uncertainty once more data is added to the likelihood function. Note that the addition of dynamic LV pressure and volume in 4 had relatively minimal effects on uncertainty when compared to only including systolic and diastolic values with 3. We recast these results into pressure-volume loops in Fig 10 and provide the 600 realizations in addition to the agglomerated average and the true model simulations. Even in the absence of atrial pressure or volume data, additional volume measurements in both the LV and RV reduce the uncertainty of atrial dynamics. LV data reduces the uncertainty substantially in 3 and 4 as shown previously in Fig 9.

Fig 9

Output uncertainty in RV and LV pressures and volumes for each experimental design.

The average model response (red) as well as ± one standard deviation (Std., gray) are provided along with the data (black circles) for each experimental design, corresponding to each column. In the first design, 1, only RV pressure is used in the likelihood, hence the uncertainty in RV volume and LV forecasts are substantially larger than that of the RV pressure. As more data is included, uncertainty in model forecasts is reduced. Note that differences between f3 and f4 are less pronounced.

Fig 10

Output uncertainty in cardiac pressure-volume loops.

Realizations in forecasts of chamber pressure-volume loops in the LA (first row), LV (second row), RA (third row), and RV (fourth row). The simplest design (1) has the largest uncertainty in simulated pressure-volume loops, except for RV pressure, which is accounted for in the likelihood. Subsequent experimental designs substantially reduce uncertainty bounds in the RV and RA (2) and eventually in the LV and LA (3 and 4).

In addition to outputs that are linked to the collected data, we investigate the uncertainty in ventricular wall strain and outcomes typically quantified during in-vivo PH studies. Engineering strain for the LV, RV, and S walls are provided in Fig 11. Strains are bounded between 5% and -20%, and there was a reduction in uncertainty when additional data was included in the likelihood. Septal strain has only a minor reduction in uncertainty for the first three designs, yet using 4 for parameter inference reduces septal strain uncertainty significantly. Moreover, using this final experimental design constrains S wall strain to have a similar shape to that of the LV and RV. Lastly, we quantify changes in mean pulmonary artery pressure, RV stroke volume, arterial and end-systolic ventricular elastance, and ventricular-vascular coupling for the different experimental designs. Histograms showing the frequency of these variables using the 600 forward samples are shown in Fig 12. Mean pulmonary artery pressure and arterial elastance have a comparable histogram width for all four experimental designs. In contrast, RV stroke volume, RV end systolic elastance, and ventricular vascular coupling have a larger variance in designs 1 and 3, which is reduced in designs 2 and 4.

Fig 11

Forecast uncertainty in LV, RV, and S wall strain.

Realizations in the LV, RV, and S engineering strain, along with the mean and ± standard derivation (Std), obtained from the posterior distributions. For designs only including RV dynamics (1 and 2), S engineering strain has a large uncertainty in the direction of strain (i.e., leftward or rightward). Designs including LV data (3 and 4) reduce the range of S strains, with the design 4 ensuring that S strain is in the same direction as the LV. LV and RV strain have substantially less uncertainty than that of S, which shrinks with more informative designs.

Fig 12

Simulated output quantities that are typically recorded when studying PH.

Histogram plots of outputs typically recorded during in-vivo studies of PH progression are generated using the same 600 samples from the posterior that were used in Figs 9, 10 and 11. These include mean pulmonary artery pressure (), RV stroke volume (SV, defined as difference between maximum and minimum RV volumes), pulmonary arterial elastance (Ea, defined the difference between and mean LA pressure divided by SV), RV end systolic elastance (Ees, defined as the end systolic ratio of RV pressure and RV volume), and ventricular-vascular coupling (VVC, defined as the ratio Ees/Ea). Differences in the experimental design had little effect on . As expected, SV was more accurately captured with additional RV volume data. The wide variability in values of Ea, Ees, and VVC using the design 1 is remedied once additional volume data is included in the design. Note that output values of VVC are made substantially more precise with additional LV data in 3 and 4.

Discussion

The present study investigates parameter identifiability for a multiscale model of cardiovascular dynamics. This work examines four different in-vivo experimental designs using in-silico modeling, and subsequently compares the reduction in parameter and output uncertainty under these different designs. In-vivo experimental designs are typically determined before using in-silico methods to analyze the data [9,30]; however, some studies have considered using the latter to plan optimal designs a-priori [45].

Sensitivity analyses

Sensitivity analyses are commonly used to reduce parameter sets to a smaller, more influential group [1,5,18]. We use these techniques to reduce the original set of 49 parameters in the model to a set of 13 influential parameters. These 13 include those attributed to the TriSeg geometry, those describing the timing, duration, and active force of sarcomere shortening, and parameters describing the systemic and pulmonary vasculature. Similar to our analysis, the study by van Osta et al. [5] used Morris screening and concluded that LV, RV, and S geometry parameters were most influential on simulations of chamber strain. While chamber strain was not considered in our experimental design, our results suggest that these same parameters are influential on ventricular pressure and volume simulations. Similar to our results, vas Osta et al. found a single parameter from the left atrium was influential [5].

We consider pressure and volume in the RV and LV as our outputs of interest, contributing to the addition of three influential circulatory parameters (R, C, and R). These three parameters were also influential in the analysis by Harrod et al. [1], who investigated PH due to LV diastolic dysfunction. The four experimental designs considered in this work focus on PH and RV function, hence more pulmonary parameters are influential than systemic. An explanation for the importance of R on RV forecasts is linked to the simplicity of the model. The total stressed volume throughout the model is held constant, hence changes in resistance or compliance will alter both pressure and volume distributions. Thus, R can have system wide effects (e.g., on RV pressure), whereas in-vivo there are mechanisms, such as the baroreflex, that can regulate system level changes in blood volume due to resistance and compliance changes. A majority of the influential parameters identified here are common in 0D models [1,18,30] and models incorporating the TriSeg framework [4,5,21], making the present analysis pertinent to future modeling studies utilizing either of these approaches.

The Morris screening methodology traditionally uses the average EE as a measure of parameter influence (5,12,35,36). While this captures which parameters are on average most influential, there may be circumstances where a parameter is highly influential in a small volume of parameter space and may require additional analyses using the maximum or median elementary effect. Previous studies have considered using the average EE for parameter fixing [5,32], yet a consistent method for parameter fixing and subset selection is warranted.

The sensitivity-based Fisher information matrix provides insight about local parameter interdependence as well as quadratic approximations of parameter confidence intervals. This method helped reduce the set of 17 influential parameters to a set of 13 locally identifiable parameters, as has been done in previous work [35]. These asymptotic analyses work well when models behave linearly within a neighborhood of the parameter value, but, as shown here with profile likelihood and MCMC, can fail in detecting practical identifiability issues.

Profile-likelihood analyses

Though local and global sensitivity analyses can identify influential parameters, they do not guarantee that parameters are identifiable [16]. The local sensitivity analysis did not reveal identifiability issues, yet the profile likelihood analysis illustrates practically non-identifiable parameters using the less detailed experimental designs. This confounding result is documented in the review by Wieland et al. [16], suggesting again that profile likelihood analyses are superior in deducing practical identifiability for nonlinear models. To the authors’ knowledge, the work by Pironet et al. [17] is the only other cardiovascular modeling study to consider this methodology. Their study [17] integrated static pressure and volume data over multiple cycles, concluding that several parameters, including total stressed volume and vena cava compliance, were practically non-identifiable. Moreover, Pironet et al. [17] reduced their initial parameter subset using sensitivity methods, but ultimately found more practically non-identifiable parameters using profile likelihood analysis. The results in Fig 5 show that six of the parameters were identifiable across all four experimental designs. Of these, five describe the structure and function of the RV and pulmonary circuit, and the last describes systemic artery resistance. Interestingly, it appears that R is practically non-identifiable using 1, but should become identifiable for larger parameter bounds. This parameter describes the state of the pulmonary vasculature, highlighting the need for additional data in the experimental design to identify its value. Parameters describing the chamber wall volumes were consistently difficult to infer, especially in the LV and S when no LV data was available. This again suggests that a true understanding of heart function requires sufficient data from both ventricles.

This study is the first to compare multiple experimental designs using a multiscale model with biventricular interaction. As expected, increasing the amount of data available reduced the confidence interval width, i.e., more data decreases the uncertainty in the estimates. In contrast to prior studies utilizing the profile likelihood [15], our results show deviations in the likelihood values for small changes in parameters. We accredit this non-smoothness to possible incompatibilities in the DAE system, which can frequently occur with this model [5,19], returning large values in the residual sum of squares. However, we are primarily interested in using the profile likelihood method to test for whether the confidence intervals have finite bounds; hence, smooth profile likelihoods are not necessary to determine if the parameter subsets are identifiable. Overall, the most complete experimental design, f4, led to the tightest confidence intervals and reduced numerical instabilities, though perturbations in V, V, and V still cause some sharp jumps in the likelihood.

Markov chain Monte Carlo

MCMC can also assess parameter uncertainty and practical identifiability [1,18,19,39]. The posterior densities in Figs 6, 7 and 8 suggest that most of the parameters are practically identifiable in the presence of measurement noise. The wall volumes, V, have wider posteriors across the first two experimental designs, but tend to shrink with additional LV data, consistent with the profile likelihood results. V has a nearly uniform posterior when only using RV pressure data (1), suggesting practical identifiability issues. This is expected, as this parameter has its largest effects on LV dynamics, which are only present in designs 3, and 4. The active stress parameter, σ, is not practically identifiable with measurement noise when using 1, and has a long posterior tail. This parameter shows noticeable changes in mixing properties when using the systolic and diastolic LV outputs in 3, and may be due to sampling in higher rejection regions to obtain appropriate LV values. All twelve pairwise plots in S2 Text using the design 4 show somewhat strong correlations between V and σ, as well A and C. This may suggest that these parameters are not practically identifiable with measurement noise. The posteriors using 3 and 4 in Figs 6, 7 and 8 are nearly all unimodal, with the true data generating parameters located near the modes. As noted by Paun et al. [39], flat, uniform posteriors suggest that parameters are not practically identifiable, supporting our claim of improved identifiability with more detailed experimental designs. The study by Harrod et al. [1] also used MCMC to test for identifiability; however, their results show a deviation between the true value of R and the posterior distribution, whereas our results (for 2, 3, and 4) show an overlap in the true and estimated values. Discrepancies between Harrod et al. and our results are attributed to the separation of systemic resistance into an arterial and venous component, whereas our model has a single systemic vascular resistance parameter corresponding to their sum. van Osta et al. [19] constructed parameter posteriors for A and the equivalent of our σ and τ using MCMC. Their study also showed that repeated construction of the posteriors from different initial guesses had reasonable overlap, suggesting all parameters were identifiable. Colunga et al. [18] contrasted two parameter subsets using MCMC and heart-transplant data. The non-identifiable set had posteriors with long, unbounded tails, whereas, like the results here, the identifiable set has tighter posterior distributions with finite tails. A comparison of the hemodynamic posteriors in Fig 8 reveals that both R and R have larger uncertainty when using the design 1. As noted previously in the text, pulmonary vascular resistance is a pertinent biomarker of PH progression and severity [23,46]. Our results suggest that, at a minimum, RV volumes are included in the experimental design to obtain reasonable estimates of hemodynamic parameters and better constrain posterior widths for V parameters. Interestingly, both the profile likelihood analysis and the MCMC results suggest that C is identifiable but without an improvement with more complex designs. This may be attributed to the simplicity of the pulmonary artery compartment, and may vary more if using a more complex model of the proximal pulmonary arteries [2].

By running multiple MCMC instances in parallel, we are able to construct individual PSRF values for each parameter and the MPSRF for each experimental design. Our results in Table 3 suggest that 50,000 iterations (with 10,000 used as burn-in) do not guarantee convergence of the MCMC process, as all of the MPSRF values are greater than 1.1. However, as detailed by Roy [44], MPSRF can be misleading in some instances. Nevertheless, we expect that substantially more iterations of MCMC (e.g., 500,000) will reduce PSRF and MPSRF values below 1.1. Our results still show that a majority of the posteriors overlap with increasing data availability in the experimental design, supporting the profile likelihood results.

Forecast uncertainty

Sampling from the parameter posteriors describes uncertainty in the model output. The first design, 1, provides information about RV pressure, and corresponding model simulations shown in Fig 9 have little uncertainty. In contrast, V(t), p(t), and V(t) exhibit larger uncertainty, with the mean response often deviating from the true signal. The more data-rich experimental designs lead to a better agreement between the model and the simulated data as well as a reduction in uncertainty. Interestingly, differences in uncertainty bounds between 3 and 4 are not evident in the isolated pressure and volume signals in Fig 9, yet pressure-volume loop uncertainty in the LV is reduced substantially in Fig 10. The difference in these two plots is linked to the timing of ventricular dynamics, which become more apparent when plotting pressure versus volume. The reduction in uncertainty when the design 3 is used suggests that including static systolic and diastolic measures of LV function are sufficient for model calibration and are necessary to reduce output uncertainty. This experimental design was utilized by Philip et al. [23] in a mouse model of PH due to left heart failure. Their results highlighted that impaired LV function can ultimately raise pulmonary vascular resistance and contribute to RV dysfunction. Assessing the LV via echocardiography is easier than the RV due to anatomic shape and location [47], hence adding this assessment to dynamic RV pressure-volume loop protocols is reasonable and provides insight into LV impairment during PH [11]. Recent studies have also found significant changes in both left and right atrial function in heart failure and PH [23,48,49]. We found only one atrial parameter, A, was influential and identifiable on RV and LV outputs. Allowing this parameter to vary explains the greater variability in left atrial pressure-volume loops than the corresponding right atrial simulations in the first three designs. However, it seems that dynamic data in the LV reduces the variability in left atrial forecasts, suggesting that 4 is the most optimal design for studying left atrial function in the absence of left atrial data. We did not consider atrial data in our possible designs, yet future work may reveal its significance in understanding disease progression, especially PH due to left heart failure [1,23].

The TriSeg model is an efficient simulator of biventricular interaction. Prior work has used this model to quantify changes in biventricular interaction under diseases such as PH [21,50], arrhythmogenic cardiomyopathy [19], and mechanical desynchrony [51]. Our results in Fig 11 show that the uncertainty in LV, RV, and S wall strain tend to decrease with more informed experimental designs. Though the model employed van Osta et al. [19] has fundamental differences from our model, both have comparable uncertainty in wall strains. Their study calibrated model predictions to measurements of wall strain by echocardiography, yet our work shows that calibration to pressure and volume data is sufficient in reducing simulated wall strain uncertainty. Strain forecasts also elucidate the state of LV-RV interaction, which is compromised in the presence of PH [11].

We use the model to simulate other hemodynamic quantities typically recorded in PH studies [23]. The distribution of simulated mean pulmonary arterial pressure in Fig 11 are similar in width across the experimental designs. Both R and C play a role in this output, yet Fig 7 shows that R has a noticeably smaller posterior when informed by 4. Though R will ultimately dictate the pressure magnitude, the unchanged posterior in C suggests that this parameter is largely attributed to mean pulmonary artery pressure. The study by Colunga et al. [18] found that including R and C in parameter inference led to close agreement between model predictions of pulmonary artery pressure and measured data. The uncertainty in mean pulmonary artery pressure described by Harrod et al. [1] are similar to our results as well. As expected, forecasts of RV stroke volume and pulmonary artery elastance (defined as the difference between mean pulmonary artery pressure and mean left atrial pressure divided by the RV stroke volume) have small variability with any designs including RV volume, i.e., 2, 3, and 4. Hence, the relatively wide probability densities for RV end-systolic elastance and RV ventricular-vascular coupling are directly tied to uncertain model predictions of RV volume. A zoom of the model forecasts shown in Fig 12 shows that additional volume constraints narrow the output uncertainty in these indices. All five indices examined here can be indicative of PH progression and RV function and suggest that RV pressure alone is not informative enough to constrain the model forecasts. Therefore, future experiments into PH and RV function should strive to have both RV pressure and volume data collected, along with static or dynamic measures of LV function.

Comparison between methods

Our results show that local and global sensitivity analyses provide insight into which parameters are influential. These two methods were the least computationally intensive; Morris’s screening with 39 variable parameters and 100 trajectories took approximately 8.7 hours, while the local sensitivity with respect to the 17 remaining parameters required 4.5 minutes of computation time. These methods only reveal whether parameters are influential and, in the case of local sensitivity, practically identifiable in the asymptotic sense via the Fisher information matrix. These methods do not guarantee that parameters are truly practically identifiable, which is where profile likelihood and MCMC analyses can be useful. However, these latter two methods are computationally expensive. Profile likelihood requires profiling a single parameter over a sufficient range with gradient based optimizations and MCMC requires numerous samples to construct the posterior. Here, profile likelihood analyses took between 30 and 50 hours depending on the experimental design, and MCMC required 110–220 hours. Neither method can be run independently in parallel, whereas sensitivity methods can be run in parallel. Nevertheless, profile likelihood analyses and MCMC uncover model and design features that cannot be identified through sensitivity analyses.

Limitations

There are several limitations in this study. The TriSeg model has been utilized by several authors to understand biventricular interaction [5,8,50]. However, this model is less detailed in handling the complex interactions between the ventricles, especially in comparison to higher fidelity finite element models. Moreover, we use diodes to represent the heart valves, which will not capture more complex dynamics seen in the tricuspid and pulmonary valve during PH [48]. A more physiological valve model could encourage echocardiographic velocity data into the experimental design. We generate synthetic data from our mathematical model to test for identifiability, hence our noise model correctly matches the true added noise. When using physiological data, this may not hold true, and may require additional components to the statistical model (e.g., model discrepancy [39]). In addition, measurement uncertainty is surely different between pressure-volume catheters and ultrasound probes and should be accounted for when using true in-vivo data for model calibration. The profile likelihood results presented here exhibit non-smoothness, whereas prior studies [17,52] typically show smooth profiles. This could be obtained by considering more sophisticated parameter mesh refinement. Our system of DAEs is stiff and can lead to model failure if parameters are not compatible. This may be overcome with more detailed information about the TriSeg geometric parameters’ covariance, which could be used to construct a non-independent prior for sampling their values during sensitivity analyses and MCMC. The posterior densities across the 12 instances of MCMC revealed that 50,000 iterations are not sufficient by PSRF and MPSRF criteria. More informative designs promoted posterior modes closer to the true parameters, yet MCMC results must be interpreted carefully if stopping criteria are not satisfied. Research into efficient MCMC and robust stopping criteria, especially in the presence of high dimensional parameter vectors, is warranted.

We consider four experimental designs that expose the coupled mechanics of the LV and RV, yet other designs could provide more insight into RV function and model calibration. A more encompassing analysis of experimental designs including additional combinations of MRI, echocardiogram, and catheter measurements in the heart chambers and vasculature is warranted. Studies using only static data will require more investigations into which parameters are most influential during systole or diastole. Our analysis is applied to data simulated for a normotensive mouse as opposed to simulating PH data. However, we believe the present analysis will be consistent even when parameters are adjusted to the PH range. This also applies to parameters dictating atrial function; our nominal model simulations do not capture the biphasic flow patterns seen in-vivo in the left and right atrium and could be included in the experimental designs in future studies. We did not consider uncertainties in volume distributions throughout the vasculature, which should be investigated further. Lastly, several parameters that require measurements at the microscale (e.g., reference sarcomere length) were fixed for our analyses. Future studies collecting data across spatial scales would require including these parameters in the above analyses and may reveal new influential parameters in the system.

Conclusion

The present study investigates parameter identifiability of a cardiovascular model with biventricular interaction, specifically calibrated for mouse hemodynamics. In summary, this study has found that:

Morris screening and local sensitivity analysis can identify influential parameters, but does not guarantee that parameters are practically identifiable;

Profile likelihood and MCMC can be utilized to identify benefits in experimental designs and deduce practical identifiability;

Model parameters describing biventricular interaction and RV function are best informed with pressure and volume data from both ventricles; and

Uncertainty in model forecasts, including cardiac pressure-volume loops and ventricular wall strain, can be substantially reduced when data from both ventricles are included.

The present analyses are conducted on model outputs corresponding to four experimental designs used to study PH and RV failure in-vivo. Profile likelihood analysis shows that model parameters are not uniquely identifiable when only RV pressure data is available, and that more informed designs are necessary to recapture the true parameter values. Our study also shows that sensitivity-based methods do not guarantee practically identifiable parameter subsets, hence profile likelihood analysis should be employed. We conclude that future, synergistic studies using both in-vivo and in-silico methods should incorporate functional LV data to improve model forecasts of cardiac function and biventricular dynamics. We hypothesize that this will be especially important when studying the progression of RV failure due to PH.

Citation diversity statement

In agreement with the editorial from the Biomedical Engineering Society (BMES) [53] on biases in citation practices, we have performed an analysis of the gender and race of our bibliography. This was done manually, though automatic probabilistic tools exist [54]. We recognize existing race and gender biases in citation practices and promote the use of diversity statements like this for encouraging fair gender and racial author inclusion.Our references contain 9.25% woman(first)/woman(last), 14.8% man/woman, 16.7% woman/man, and 59.3% man/man. This binary gender categorization is limited in that it cannot account for intersex, non-binary, or transgender people. In addition, our references contain 3.70% author of color (first)/author of color(last), 5.55% white author/author of color, 25.9% author of color/white author, and 64.8% white author/white author. Our approach to gender and race categorization is limited in that gender and race are assigned by us based on publicly available information and online media. We look forward to future databases that would allow all authors to self-identify race and gender in appropriately anonymized and searchable fashion and new research that enables and supports equitable practices in science.

Summary of model parameterization and initial conditions, full results from Morris screening, additional MCMC results, and list of model parameters.

(DOCX)

Click here for additional data file.

Pairwise plots from all four experimental designs and all 12 MCMC instances.

(DOCX)

Click here for additional data file.

Comparison plot of the modified average elementary effect, μ*, against elementary effect variance, s2, for each of the output quantities of interest.

Results from the Morris screening show a trend of higher mean effects and higher variance in elementary effects, which indicate an influential parameter.

(EPS)

Click here for additional data file.

Relative error in the median posterior values compared to the true, data generating parameters.

Black squares represent the true parameter values while black circles indicate the mean value across the median posterior value for the 12 instances of MCMC.

(EPS)

Click here for additional data file.