An in-silico analysis of experimental designs to study right ventricular function and pulmonary hypertension pulmonary in mice: a model-based analysis of Analysis of Right Ventricular Mechanoenergetics.

14 In-vivo studies of pulmonary hypertension (PH) have provided key insight into the progression of the disease and right ventricular (RV) dysfunction. Additional in-silico experiments using 16 multiscale computational models have provided further details into biventricular mechanics and 17 hemodynamic function in the presence of PH, yet few have assessed whether model parameters 18 are identifiable prior to data collection. Moreover, none have used modeling to devise synergistic 19 experimental designs. To address this knowledge gap, we conduct an identifiability analysis of a 20 multiscale cardiovascular model across four simulated experimental designs. We determine a set 21 of parameters using a combination of Morris screening and local sensitivity analysis, and test for 22 identifiability using profile likelihood based confidence intervals. We employ Markov chain 23 Monte Carlo (MCMC) techniques to quantify parameter and model forecast uncertainty in the 24 presence of noise corrupted data. Our results show that model calibration to only RV pressure 25 suffers from identifiability issues and suffers from large forecast uncertainty in output space. In 26 contrast, parameter and model forecast uncertainty is substantially reduced once additional left 27 ventricular (LV) pressure and volume data is included. A comparison between single point 28 systolic and diastolic LV data and continuous, time-dependent LV pressure volume data reveals 29 that even basic, functional data from the LV remedies identifiability issues and provides 30 substantial insight into biventricular interactions. years after disease onset. Many researchers couple computational models with in-vivo experimental models of PH, yet few ever assess what data might be necessary or sufficient for parameter inference prior to designing their experiments. Here, we considered a multiscale computational model including sarcomere dynamics, biventricular interactions, and vascular hemodynamics, and assessed whether parameters could be inferred accurately given limited cardiac data. We utilized sensitivity analyses, profile likelihood confidence intervals, and MCMC to quantify parameter influence and uncertainty. We observed that RV pressure alone is not sufficient to infer the influential parameters in the model, whereas combined pressure and volume data in both the RV and LV reduced uncertainty in model parameters and in model forecasts. We conclude that synergistic PH studies utilizing computational modeling include these data to reduce issues with parameter identifiability and minimize uncertainty.

6 105 monitoring animal models of PH and focus on the RV (23-25), and an additional design that 106 utilizes dynamic pressure-volume data in both the LV and RV (26). We generate both noise-free 107 and noisy data from the model to test for parameter identifiability and analyze the output 108 uncertainty in model simulations and several biomarkers of PH progression.   (2) 8 146 (mm 2 ), and radius of mid-wall curvature (mm -1 ) (8). Once has been calculated and the 147 corresponding is obtained from the sarcomere model, the mid-wall tension can be calculated 148 as = 2 ( 1 + 2 3 + 4 5 ) . ( . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 25, 2022. ; https://doi.org/10.1101/2022.03.22.485347 doi: bioRxiv preprint 158 where , ( l) is the unstressed volume (see S1 text), ( l KPa -1 ) is the vascular compliance, 159 and (KPa s l -1 ) is the vascular resistance between compartments. Finally, we model the two 160 atrioventricular valves (mitral and tricuspid), the two semilunar valves (aortic and pulmonic), 161 and the resistor between the systemic veins and right atrium as diodes 162 The hemodynamics model consists of eight differential equations for ( ), eight resistance 163 parameters, and four compliance parameters.
164 Summary. The multiscale model consists of 18 differential equations (describing , Γ , and ), 165 two algebraic constraints (equation (4)), and a total of 49 parameters. Due to the algebraic 166 constraints, the model constitutes a system of differential algebraic equations (DAEs) and is 167 solved using the variable-step, variable-order ode15s solver available in MATLAB (Mathworks; 168 Nantick, MA).

175
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

187
The local sensitivity of a model output with respect to a parameter, , is approximated 188 by the centered difference 189 where ( ; ) is the quantity of interest from the model, Δ is the step change in parameter 190 value, and is the -th unit vector. For time-dependent outputs, we consider the 2-norm of the 191 model output, i.e. , = | , | 2 2 . We account for differences in parameter magnitude by 192 computing the log-scaled parameter sensitivity, / log( ) (31,35).

193
The Morris' screening approach computes the "elementary effects" . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made 194 where ( ) is the parameter step size describing the "levels" of effects. Choosing to be even 195 provides a more symmetric sampling distribution (33), hence we choose = 60 giving ≈ 0.51.
196 Note that EE , is a coarser approximation of model sensitivity than , , but is quantified over a 197 larger parameter space. We scale parameters from their original value to the interval [0,1] as 198 done previously (34), and utilize the algorithm provided by Smith (12) to construct our sampling 199 methodology. The indices from the Morris method are 200 Here, , is the average of EE , , * , is an improved metric for average model sensitivity (34), 201 and 2 , is the variance of EE , . We use the combined index, , = * , 2 + 2 , , to measure a 202 parameter's influence (36).

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

206
207 Practical parameter identifiability 208 In this work, we assess parameter identifiability using three techniques. The first is 209 through the local and global sensitivity metrics discussed above. Next, we consider the profile 210 likelihood, which provides information about whether each is identifiable from a given set of 211 data. Lastly, we use MCMC methods for Bayesian inference, and utilize the marginal posterior 212 distributions to assess parameter identifiability.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made 213 Sensitivity based identifiability. Parameters that have little effect on the model output are 214 considered non-identifiable, since they do not affect the quantity of interest (12), and should be 215 fixed before conducting inference. We employ a two-part parameter fixing methodology using 216 the results from Morris screening and local sensitivity analysis.

217
A parameter is deemed non-influential for all outputs if its index , is less than the 218 average for all parameters = 1,2,…, 219 where is one of the model outputs (5,32). Parameters that are less than this threshold for all 220 outputs are considered non-influential for inference and are fixed.

221
After using the Morris screening approach, the subset is analyzed by conducting a local 222 sensitivity analysis around the nominal parameter values. The Fisher information matrix, = ⊤

223
, must be non-singular for gradient based parameter estimation, hence its utility in parameter 224 identifiability (37). If is invertible but has a large condition number, then some of the 225 sensitivities are nearly linearly dependent and the subset requires further reduction. We use an 226 eigenvalue-eigenvector analysis method to determine which parameters cause the ill-227 conditioning of (14,38), and fix these parameters at their nominal value.
228 Profile likelihood. The most common and robust technique for assessing practical identifiability 229 is the profile likelihood (13,15). This technique increments a fixed parameter, , while 230 minimizing the negative log-likelihood for all other parameters in the subset, i.e.
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made 231 where is the -th data source, is the corresponding model output, LL( | ) is the log-232 likelihood, 2 is the noise variance for the data source, and is the number of data points. The 233 corresponding profile likelihood confidence intervals for are (13) 234 Each ( ) is constructed around the optimal estimate, * , and depends on the inverse 235 cumulative distribution function of the chi-squared distribution, ( 2 1 , ), with one-degree of 236 freedom and confidence level (13). If ( ) is completely flat (e.g., ( ) is infinite), then 237 is deemed structurally non-identifiable and cannot be uniquely determined due to model 238 structure. If only one side of ( ) is flat, then is considered practically non-identifiable, and 239 could become identifiable if more data was available for inference (16).
240 Bayesian inference. Bayesian parameter inference using MCMC is more computationally 241 expensive than gradient based inference, but provides detailed insight into parameter 242 relationships and avoids local minima in the likelihood (39-41). We use the DRAM algorithm 243 (42), which is described in depth elsewhere (31,43). In short, the goal of MCMC is to 244 approximate the posterior distribution

254
We utilize DRAM on a set of noisy data, generated by the model and corrupted with 255 noise. To ensure adequate parameter space coverage, we first initialize a gradient based 256 optimization from twelve randomized initial parameter sets and minimize the residual sum of 257 squared errors for the given experimental conditions (defined in the next section). The optimal 258 parameter vector, , is used as a starting value for DRAM and the Hessian matrix obtained 259 from the optimization is used as the initial covariance matrix to preserve possible sampling 260 asymmetry (12). We implement this using the freely available DRAM package developed Haario . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made

281
We perform all sensitivity and identifiability analyses with respect to the pressure and (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made We ran the Morris screening algorithm using 100 randomized initializations. Fig 3 shows 306 the parameter ranking , using the mean effect * and corresponding variance 2 for the RV 307 and LV pressures and volumes. See S1 text for individual results from the Morris screening 308 analysis as well as parameter bounds for sampling. Sensitivity results were analyzed by 309 comparing the parameter ranking , to the mean effect for each ventricular pressure and 310 volume. All four compliances were consistently ranked within the most influential parameters, 311 while other parameters describing cardiac chamber dynamics (e.g., , ) varied with the 312 output. We fixed parameters that were less influential than for all four outputs (i.e., pressure . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made , , . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  Table 2, was used in the profile likelihood and MCMC analysis.

334
Profile likelihood based confidence intervals are constructed using the noise-free, model 335 generated data. We construct the confidence intervals ± 50% away from the true parameter 336 value, with the confidence level cutoff for each design calculated using equation (14)   Large deviations in the profile likelihood correspond to local minima and parameter sets that are 351 incompatible for the system of DAE's.

352
. CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 25, 2022. ; https://doi.org/10.1101/2022.03.22.485347 doi: bioRxiv preprint Noise corrupted data generated by the model is used in the likelihood defined in equation 354 (13). We use minimally informative priors (i.e., with a large variance) for each parameter and 355 initialize the DRAM algorithm using the optimal parameter vector and estimated covariance 356 matrix from twelve randomly selected initial guesses. MCMC is run for 50,000 iterations, with 357 the initial 10,000 being left out as a "burn-in" period. We separate the results from MCMC into 358 three groups: parameters representing the heart chambers' geometry (Fig 5), parameters within 359 the sarcomere model (Fig 6), and hemodynamic parameters in the circulatory model (Fig 7).

375
The model parameters indicative of the TriSeg geometry (wall volume, , and reference mid-

423
The simplest design ( 1 ) has the largest uncertainty in simulated PV loops, except for RV 424 pressure, which is accounted for in the likelihood. Subsequent experimental designs substantially 425 reduce uncertainty bounds in the RV and RA ( 2 ) and eventually in the LV and LA ( 3 and 4 ).

427
In addition to outputs that are linked to the collected data, we investigate the uncertainty 428 in ventricular wall strain and outcomes typically quantified during in-vivo PH studies.

451
Histogram plots of outputs typically recorded during in-vivo studies of PH progression are 452 generated using the same 600 samples from the posterior that were used in Fig. 8, 9, and 10. in-silico methods to analyze the data (9,30); however, some studies have considered using the 470 latter to design optimal designs a-priori (44).

471
Sensitivity analyses are commonly used to reduce parameter sets to a smaller, more 477 They concluded that LV, RV, and S geometry parameters were most influential on simulations of 478 chamber strain. While chamber strain was not considered as one of our output quantities of 479 interest, our results suggest that these same parameters are influential on ventricular pressure and 480 volume simulations. Their results also showed that, similar to our results, a single parameter 481 from the left atrium was influential (5). We consider pressure and volume in the RV and LV as . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made 505 structure and function of the RV and pulmonary circuit, and the last describes systemic artery 506 resistance. Interestingly, it appears that pulmonary vascular resistance, , is practically non-507 identifiable using 1 , but should become identifiable for larger parameter bounds. This parameter 508 is often used to describe PH and the state of the pulmonary vasculature, highlighting the need for 509 additional data in the experimental design. Parameters describing the chamber wall volumes 510 were consistently difficult to infer, especially in the LV and S when no LV data was available.
511 This again suggests that a true understanding of heart function requires sufficient data on both 512 ventricles.

513
This study is the first to compare multiple experimental designs using a multiscale model

525
MCMC is a common tool for Bayesian inference, and can also assess parameter 526 uncertainty and identifiability (1,18,19,39). The posterior densities in Fig 5, Fig 6,  530 Specifically, , has a nearly uniform posterior when only using RV pressure data ( 1 ), 531 suggesting practical identifiability issues. This is expected, as this parameter has its largest 532 effects on LV dynamics, which are only present in designs 3 , and 4 . The active stress 533 parameter, , , is not identifiable using 1 , and has a long tail towards larger values. This 534 parameter shows noticeable changes in mixing properties when using the systolic and diastolic 535 LV outputs in 3 , and may be due to sampling in higher rejection regions to obtain appropriate 536 LV outputs. All twelve pairwise plots in S2 text using the design 4 show somewhat strong 537 correlations between , and , , as well , , and . Contrasted with the profile 538 likelihood results, this may suggest that these parameters are not practically identifiable given 539 measurement noise. The posteriors using 3 and 4 in Fig 5, Fig 6,  580 Assessing the LV via echocardiography is easier than the RV due to anatomic shape and location 581 (46), hence adding this assessment to dynamic RV PV loop protocols is reasonable and provides 582 insight into LV impairment during PH (11). Recent studies have also found significant changes 583 in both left and right atrial function in heart failure and PH (23,47,48). We found only one atrial 584 parameter, , , , was influential and identifiable on RV and LV outputs. Allowing this 585 parameter to vary explains the greater variability in left atrial PV loops than the corresponding 586 right atrial simulations in the first three designs. However, it seems that dynamic data in the LV 587 reduces the variability in left atrial forecasts, suggesting that 4 is the most optimal design for 588 studying left atrial function in the absence of left atrial data. We did not consider atrial data in 589 our possible designs, yet future work may reveal its significance in understanding disease 590 progression, especially PH due to left heart failure (1,23).

591
The

682
The sarcomere length is calculated as a function of the myocardial strain

689
The mechanical activation of the sarcomere is heuristically modeled as separate "rise" 690 and "decay" functions. The former is given by . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made  where (KPa) is a parameter that scales the active force contribution.

702
For the passive stress, we consider a new, passive reference length , , ( m), giving 703 the passive strain . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made 715 which is updated at each time point. Hence, we solve for using in equation (5)   CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 25, 2022. ; https://doi.org/10.1101/2022.03.22.485347 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 25, 2022. ; https://doi.org/10.1101/2022.03.22.485347 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made The copyright holder for this preprint this version posted March 25, 2022. ; https://doi.org/10.1101/2022.03.22.485347 doi: bioRxiv preprint . CC-BY 4.0 International license available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made