Convergence, sampling and total order estimator effects on parameter orthogonality in global sensitivity analysis

Dynamical system models typically involve numerous input parameters whose “effects” and orthogonality need to be quantified through sensitivity analysis, to identify inputs contributing the greatest uncertainty. Whilst prior art has compared total-order estimators’ role in recovering “true” effects, assessing their ability to recover robust parameter orthogonality for use in identifiability metrics has not been investigated. In this paper, we perform: (i) an assessment using a different class of numerical models representing the cardiovascular system, (ii) a wider evaluation of sampling methodologies and their interactions with estimators, (iii) an investigation of the consequences of permuting estimators and sampling methodologies on input parameter orthogonality, (iv) a study of sample convergence through resampling, and (v) an assessment of whether positive outcomes are sustained when model input dimensionality increases. Our results indicate that Jansen or Janon estimators display efficient convergence with minimum uncertainty when coupled with Sobol and the lattice rule sampling methods, making them prime choices for calculating parameter orthogonality and influence. This study reveals that global sensitivity analysis is convergence driven. Unconverged indices are subject to error and therefore the true influence or orthogonality of the input parameters are not recovered. This investigation importantly clarifies the interactions of the estimator and the sampling methodology by reducing the associated ambiguities, defining novel practices for modelling in the life sciences.


Introduction
Parameter identifiability addresses the question of whether, and to what degree it is possible to uniquely estimate input parameters (inputs) for a given dynamical system with a set of measured (or synthetic) outputs.This problem is typically decomposed into practical identifiability, which incorporates practical estimation issues associated with real data (such as noise and bias), and structural identifiability, which considers only model structure.The latter investigation is deemed ideal and effectively assumes that all data are known at every time point and are free of bias and noise [1].Practical identifiability accounts for the role of noise and sampling frequency inter alia in hindering the ability uniquely to estimate inputs.These issues notwithstanding, the study of unique parameter estimation is very important to the complex models increasingly used in life sciences, which encompass pharmacology, epidemiology and cardiovascular applications [2][3][4].
Assuming one can identify inputs representative of the data, we arrive at model personalisation-a process of effectively calibrating a life science model using data available from an individual subject or patient.Within a clinical setting, this might involve calibrating a cardiovascular model to (patho)physiological metrics.Robust and reliable model personaisation is a key component for the development of digital twins for healthcare applications [5].
Unique model parameters are normally obtained by solving an inverse problem.One seeks the extrema of a cost function, typically the L 2 norm of some weighted difference between measured, often very noisy and self-inconsistent target experimental data and a corresponding model prediction.A model cost function interacts with a hyper-surface in the model's input parameter space and personalisation amounts to an attempt to locate the global minimum of the cost function.It is appropriate to emphasise, here, that input parameter space is multidimensional (with a dimension equalling the number of input parameters) and that the gradient of the cost function hyper-surface typically varies rapidly in some parameter directions (axes), but very slowly in others.Hence, a sufficient consideration of a model's input parameter hyper-surface is central to understanding which parameters can be recovered uniquely.
Various methods exist to calculate model identifiability (see e.g [6]); our approach is based on the method of Li et al. [7], which calculates the identifiability index of the ith model input I i as follows: where d i is the ith input's orthogonality relative to a pre-selected, existing set of input parameters (which one is seeking to expand) and E i is an ith input's effect.Here, the index I i measures the likelihood of a unique recovery of the ith model input.In our method, both effect and orthgonality of an input are calculated from the sensitivity matrix, generated with respect to model outputs [7][8][9][10].Clearly, the identifiability index depends on both the effect and the orthogonality which can be disclosed by sensitivity analysis.This prompts an investigation to find the most reliable and most robust method for calculation.
Sensitivity analysis studies how a change in a model's specific output can be attributed to different sources of uncertainty in its (likely) many inputs [11].Two types of sensitivity analysis exist: (i) Local sensitivity analysis (LSA), which examines the sensitivity of the model inputs at one specific model operating point in the input parameters' space; (ii) Global sensitivity analysis (GSA), which determines sensitivities at multiple points throughout input space, before finding a statistic measure [12].Variance-based indices are commonly recognised as the pre-eminent statistic of GSA, due to their model-independent nature, and their ability to account for interactions between model inputs and the ease of interpretation [13,14].For studies on practical identifiability of inputs, the metrics on the total order sensitivity matrix are calculated following Eq (1), which evaluates the overall contribution of an input and its interaction with other inputs to a specific output whilst considering orthogonality.We defer further detailed discussion to Section 2.3.
Calculations of inputs' effects and orthogonality are an important area of research for model personalisation.Here we ask three universal and interrelated questions: 1. What is the most reliable estimator for the underlying sensitivity indices to be computed on?
2. What is the optimal sampling methodology, in relation to (i) above, for which one explores the complex input parameter space?
3. How do the choices of estimator and sampling methodology impact the index convergence?
Previous work [15][16][17] mostly concentrates on efficient and accurate computation of the total order matrix and the evaluation of different estimators' abilities to reveal the "true" effects of inputs, given their interactions.Recently, Puy et al. [18] reported their examination of several total order estimators-essentially a sensitivity analysis of a sensitivity analysis [19].These authors varied the sampling method, between Monte Carlo and quasi-Monte Carlo, their analytic (note) test model, the dimensionality of input parameter hyperspace, the distribution of input parameters, and the number of model runs.The work provides a comprehensive and systematic assessment of the properties of different estimators.
We structure the paper as follows: Section 2 introduces the sampling methodologies and the total order estimators used for our investigation, as well as the nonlinear systems they are applied to; Section 3 presents our findings and section 4 declares and discusses best practice to offset interactions between sampling methodologies, estimators and model dimensionality when considering the orthogonality of inputs.

Backgound
Within the field of parameter identifiability, it is accepted [20] that there is interaction between the methodology of sampling from the model's input parameter hyperspace and the estimators used to extract values of, e.g.Sobol indices [16].It follows that other sensitivity derived metrics, such as the parameter orthogonality [7], are also affected by this same interaction.In this work, therefore, we extend the consideration of both sensitivity and orthogonality in tandem.We choose to undertake this assessment based upon a very important class of cardio-vascular models, which are notoriously stiff and non-linear, in contradistinction to the simple, traditional models usually used for such purposes [18,19,21], which possess analytic solutions.Because there is no analytical solution to the models considered in this work, the reference (or true) values of the sensitivity and orthogonality indices must therefore be determined using high levels of (re)sampling, based upon accepted processes of resampling and bootstrapping [21].
Our dynamical system has the form: where Y ðt; yÞ represents the dynamic outputs described by the measurement functions, h , which depend on the state variables X ðt; y Þ.The state variables are parameterised by the inputs y but (for our application, at least) h does not directly depend upon y .f represents a set of square integrable functions which, together with inputs, y , determine X ðt; y Þ, the solution of equations ( 2).Accordingly, f ðX ðt; y Þ; y ÞÞ, now viewed as a function of y , is a suitable surrogate to determine the sensitivity of model inputs on outputs, through X ðt; yÞ.For inputs which range over a bounded region, we consider a functional I½ f �ðtÞ, formed from the integral of f over the unit hypercube I n = [0, 1] n , where n represents the dimensionality of the input parameter space: The effect of sampling inputs can be assessed with reference to this integral, by viewing the sampling process as an effective quadrature in which the sampled inputs define the abscissae.
The quality of the sampling can then be measured by the quality of the quadrature.The important question of how its accuracy is determined in conjunction with the sampling of the hypercubic region of input parameter space is central to a robust and reliable sensitivity analysis.

Sampling methodologies
In this section, we declare and briefly describe the input parameter space sampling methodologies that will be assessed in this work.We concern ourselves with two popular Monte Carlo (MC) sampling methodologies: Uniform (U) and Latin Hypercube (LH), and three Quasi-Monte Carlo (QMC) sampling methods: Golden Ratio (GR), Lattice Rule (LR) and Sobol Sequence (SS).

2.2.1
Monte-Carlo sampling methods.Uniform sampling.The simplest sampling approach from literature is uniform sampling [22].Input parameters y are regarded as uniformly distributed random variables, within the hypercube I n such that: where E is the expectation operator, y is a parameter vector of length n.This is deemed a crude approximation with poor convergence rates [23].Latin hypercube sampling.The efficiency of MC methods is determined by the properties of the random samples.A priority for researchers is to develop strategies which ensure points are placed more uniformly, within I n .One response is to use LH [24], which is a very common methodology in life sciences [25][26][27][28].Its main objective is to reduce the variance associated with evaluating Eq (3).One decomposes the space of inputs into N-dimensional squares, to ensure the space is sampled as uniformly as possible.Let f� i j g, for j = 1, . .., n, be independent random permutations of samples i = 1, . .., N, each uniformly distributed over all N! possible permutations.One then sets where U i j are independently randomly sampled points on the [0, 1] interval.It can be seen intuitively how only one sample point of the input parameters falls between iÀ 1 N and 1 N for each dimension j = 1, . .., n.Here, we use the 'standard' version of LH, but see [29] for other variations.
2.2.2 Quasi Monte-Carlo sampling methods.An improvement on the Monte-Carlo sampling methodologies is the low discrepancy sampling (LDS) methods, coupled with the QMC algorithm, as shown in [30,31].Discrepancy is a measure of the deviation of sampled points from the uniform distribution [32].Consider a number of points N R from a sequence {θ i }, for i = 1, ‥, N, in an n-dimensional rectangle R centred upon an origin 0, whose sides are parallel to the coordinate axis, which is a subset of I n : R � I n , where R is attached with a measure.A sequence has low discrepancy if the proportion of points in the sequence falling into an arbitrary set R is close to the measure of R. LDS satisfies the upper bound condition [33]: where D N is the sample discrepancy and k(n) is a particular constant depending on the sequence and size of input paraeter space.LDS is designed to place sample points as uniformly as possible mathematically, within a hypercube, instead of the statistical approach adopted in LH.The QMC approximation of the integral in Eq (3) has identical form to Eq (4).
except in this framework, y i;q is a quasi sampled parameter vector which has been generated from an LDS and the points are distributed uniformly in the unit hypercube I n .As a consequence, the sample points generated on I n have a deterministic nature.Golden ratio sampling.Golden ratio sampling is an LDS sampler in which sample points are based on the fractional part of successive integer multiples of the golden ratio.First introduced by Schretter and Kobbelt [34], using a simple incremental permutation of a generated golden ratio sequence, they demonstrated equal coverage of a two-dimensional space.Since this point experimental investigations have begun to explore the effectiveness of this sampling method [35].Let y ¼ ðy 1 ; y 2 ; . . .; y n Þ be a n-dimensional parameter vector, where each parameter θ i has a lower bound y min i and an upper bound y max i .We want to generate N samples of y .First, we define a set of n distinct points fx i g n i¼1 in R n as follows: Where x i is an alternating sequence of cosines and sines.Then, we define a set of n weights fw i g n i¼1 using the golden ratio ϕ: To generate N samples fy j * g N j¼1 , we use the following procedure: where � denotes the element-wise product, y min and y max are the vectors of lower and upper bounds, respectively.Samples are then perturbed as in section 2.2.1 to generate N independent samples.It is the structure obtained from utilising the golden ratio that provides the effective structure.In this work and specifically in Section 3, we will test GR sampling for systems with input parameter dimensions much higher than two.
Rank-1 lattice rule sampling.Another LDS is the rank-1 lattice rule, where an n-dimensional rank-1 lattice P is a set of points that contains no limit points and satisfies [36]: y 0 2 P ) y þ y 0 2 P and y À y 0 2 P; 8y : ð11Þ A general lattice is constructed by a generating matrix G 2 R n�n : where V is any integer unimodular vector.A generator matrix is not unique to a lattice P, i.e., P can be obtained from different generator matrices.A rank-1 lattice is a special case of the general lattice, which has a simple operation for point set construction, instead of directly using Eq (12).A rank-1 lattice point set can be constructed as where z 2 Z n is the generating vector and the inner product denotes the operation of taking the fractional part of the input number element-wise.One can then scale y i as above to obtain parameter samples between the specified bounds.Compared with the general lattice rule, the construction form of the rank-1 lattice already ensures the constructed points are inside the unit cube, without the need for any further checks.Sobol sequence sampling Our final sampling methodology is the well-known Sobol LDS [37].The Sobol sequence is widely considered as the optimal sequence for exploration of an input parameter space [27,[38][39][40].The Sobol low-discrepancy sequence is a quasi-random sequence based on the following equation: a i,j is the j-th digit in the binary representation of i, and v j is the j-th direction vector in n dimensions.
To ensure that the samples y i lie within specified bounds for each dimension, we can modify the equation as follows: The direction vectors v j are pre-computed and can be obtained from various sources.One common choice is to use the direction vectors provided by Sobol [37].Put simply, the Sobol LDS aims to achieve three requirements: (1) best uniformity as n ! 1 as often sampling methodologies struggle to provide adequate coverage in high dimensional spaces [41][42][43]; (2) good distribution even with small parameter sizes; (3) a very fast computational algorithm.

Sensitivity analysis and orthogonality
Given a model of the form in Eq (2) with Y (a continuous or discrete output), a variance based first order or total order effect can be calculated for a generic input factor θ i .y c i denotes the complementary set, i.e., all other model inputs excluding θ i .The first order sensitivity index can be written as: where E is the expectation operator.The inner expectation operator functions such that the mean of Y is taken over all possible values of y c i while keeping θ i fixed.The outer variance is taken over all possible values of θ i .Then utilising the known identity [44]: where Var y i ðE y c i ðYjy i ÞÞ measures the first order (additive) effects of θ i on the model outputs.
Another popular variance measure (and the concentration of this work) is total order estimators, first introduced by Homma [45]: Here S T,i measures the total effect, i.e., first and higher order effects (multiplicative interactions) of input parameter θ i .One can consider this by recognising that Var y c i ðE y i ðYjy c i ÞÞ is the first order effect of y c i , so VarðYÞ À Var y c i ðE y i ðYjy c i ÞÞ must give the contribution of all terms in the variance decomposition which do include the input θ i .
The equations can be derived through a Hoeffding-Sobol decomposition, and utilising the fact that each term is assumed to be square integrable.The detailed derivation can be found in [46,47].
Sobol indices converge slowly in general and estimators are used to accelerate the process [11,13].Here, we benchmark five sampling methodologies against four commonly chosen total order estimators: Homma and Saltelli [45], Sobol [48], Jansen [49] and Janon et al. [50].While this list is far from exhaustive, it represents a selection of total order estimators which have been practically used within the field and are not costly to execute computationally.
For the Homma & Saltelli, Sobol and Janson estimators, their mean and variance take the following form: and for the Janon estimator: In order to utilise the above formulae, we propose two independent sampling matrices are generated-A and B, with elements a ij and b ij , for i = 1, . .., N, j = 1, . .., n (where N is the number of samples and n is the total number of input parameters).We can now introduce a matrix A i B or B i A where all the rows are from A or B, except the i-th row which is extracted from B or A. These matrices are then used to compute the sensitivity indices which will be discussed below.While one can use either combination of the A and B matrices, within this work we use the couple A; A i B due to its proven efficiency in calculating sensitivity indices [16].When utilising the estimators in Table 1, f represents the differential algebraic equations (DAE) system which we solve numerically.f is then executed by supplying a parameter vector (columns from our sample matrices A and B).We then obtain a total order sensitivity matrix S T , of dimension (m × n) where m and n are the measurement and parameter dimensions, respectively.
As will be defined in the next section, we will calculate total order indices for both continuous and discrete measurements.For continuous measurements, calculating the total order index produces waveform data which demonstrate the sensitivity of each input parameter over the cardiac cycle.In order to quantify the effects continuous measurements have on the calculation of total order, we must average this sensitivity waveform.Rather than averaging across a time range (which process regions of low variance equally to those of high variance), we seek to expose differential sensitivities by examining variance-weighted averages: where TAS T,i is the time averaged total order effect of an input parameter i and Y c ðt k Þ represents the approximated continuous measurement at time step k.The division is performed component-wise such that TAS T,i is of size (m × 1) and this vector represents the time averaged sensitivity indices for an input parameter i against all outputs m.We use the measure d ij to measure orthogonality between input parameters θ i and θ j : where S T T;j represents the transposed total order sensitivity vector of input parameter j.The standard dot is used to compute S T T;j � S T;i and the Euclidean norm is utilised on each sensitivity vector ||S T,i ||.Calculating Eq (22) returns a n × n matrix where each element d ij 2 [0, 1].d ij = 0 represents total dependence between input parameters and d ij = 1 represents complete independence.We then are able to rank input parameters by calculating the mean of input parameters' orthogonality and ranking them such that the input parameter with the highest orthogonality is in position 1 and the input parameter with the lowest orthogonality is in position n, where n is the input space dimension.

Models and data
We examine two models, which are representative of the heart and systemic circulation, of varying dimensionality in order to assess their potential effects on total order estimators.The Table 1.Formulae to compute S T , where f 0 and Var represent the mean and variance of the outputs respectively, as defined in Eqs ( 19) and (20).l runs from 1 to N for the number of model samples.

Authors
Estimator S T Homma & Saltelli [45] VarðYÞ À Janon et al. [50] VarðYÞ À All model simulation code is avilable at https://github.com/H-Sax/Orthgonality-SA.The latter model is of the same form as Bjordalsbakke's, however, the input parameter space dimensionality has increased from nine to twenty, due to the addition of an atrium and other compartments.
Each compartment state is specified by its dynamic pressure P(t) (mmHg), an inlet flow Q (t) (mL/s) and a volume V(t) (mL): where i 2 {lv, sa, sv} for the 1-chamber model and lv, sa, sv represent the left ventricle, systemic artery and systemic veins respectively; i 2 {lv, la, sas, sat, sar, scp, svn} for the 2-chamber model and lv, la, sas, sat, sar, scp, svn denote the left ventricle, left atrium, systemic aortic sinus, systemic artery, systemic arterioles, systemic capillary and systemic vein.Formally, t is the continuous variable time.The input parameters for the 1 and 2 chamber models are displayed in Tables 2 and 3.
In generic form, the equations relating to the passive compartmental state variables all take the form: where the subscripts (i − 1), i, (i + 1) represent the proximal, present and distal system compartments, respectively.V s,i (mL) denotes the circulating (stressed) volume [53].C i (mL/ mmHg), R i (mmHgs/mL) and L i (mmHg s 2 /mL) denote compartmental compliance, the Ohmic resistance and compartmental inertia between compartments i and (i + 1).See Fig 1 and Tables 2 and 3.
where R val represents the resistance across the respective valve.
Let us consider the active model compartment.The dynamics of the left ventricle or left atrium are described by a time-varying compliance C(t), or reciprocal elastance, E(t) (mmHg/ mL) which determines the change in pressure for a given change in the volume [53]: where V 0 & V s (t) represent the unstressed and stressed volumes, respectively, in the left ventricle or left atrium.E(t) may be described in analytical form as follows [52]: where e(t; τ es , τ ep ) is the activation function for both the ventricle and the atrium and is parameterised by the end systolic and end pulse timing parameters τ es and τ ep respectively.The elastance function is defined over one cardiac cycle, i.e., time � t 2 ½0; t� with τ (the length of the cardiac cycle) fixed in this work to τ = 1 s.The contractility, E max , and the compliance, E min , both control the elastance extrema of the left ventricular and the left atrium.There is a discontinuity in E(t) at t = τ when the next cycle starts.
Model A is implemented directly as a system of ODEs, Model B is implemented using an acausal modelling framework which will simplify the implementation of model variations.The acausal modelling library is published as a Julia package, CirculatorySystemModels.jl [55].For further information on the mathematical construction of these systems, see [56].
As we focus on the computation of total order indices, along with varying model dimensionality, we will also vary the type of data provided to the model between continuous and discrete.For the simple nine parameter model defined in Fig 1A , we utilise: where Y c and Y d represent the continuous and discrete measurement vector.For the twenty dimensional model in Fig 1B we utilise similar measurements: The discrete measurements above are scalar quantities extracted from a continuous waveform solution.The mean measurement averages across the whole waveform which then gives the scalar measurement.
We follow the advice of Saltelli et al. [16] that the first and total order indices require k computations: where n is the dimensionality of the input parameter space, N is the number of samples taken from the input parameter space and k is the total number of model evaluations needed in order to compute the indices.It is well accepted within the literature, that utilising the low discrepancy Sobol sequence combined with the Jansen Estimator is considered best practice for calculation of total order indices [16,18].Given the high non-linearity and stiffness associated with these lumped parameter models, we first perform an assessment of convergence and uncertainty in the calculation of the total order indices.We consider the indices converged once there is no more change in the rank of the input parameters and the error associated with the indices is less than 5%.Once the indices have converged with the Jansen estimator and Sobol sampling methodology, we fix this sample size, N, for all other estimators and sampling methodologies.The uncertainty associated with the indices is calculated using re-sampling with replacement, where we set the number of bootstraps to B = 1000, as found in the literature [21,57].Once convergence has been achieved for the discrete measurements, we use this sample size N for the continuous outputs to ensure the time averaged indices derived are not subject to excessive uncertainty.
All computations are performed using Julia [58] and reproducible code is available at https://github.com/H-Sax/Orthgonality-SA.Step 1 involves defining the ODEs of the system of interest.Model A is defined employing the package DifferentialEquations.jl[59] and Model B is implemented using our acausal modelling library CirculatorySystemModels.jl [55].Step 2 generates the sample points needed to perform GSA.QuasiMonteCarlo.jl is a Julia package which provides the needed sampling algorithms.Step 3 uses GlobalSensitivity.jl [60] which is a Julia implementation of the Sobol indices with varying estimators and bootstrapping methodology.Then to analyse the sensitivity indices, we utilise self written functions.
Specifically, simulations were solved using Vern7 algorithm [61], with relative and absolute tolerances set to 10 −8 .We saved the model solution at 200 points between cycles 15 and 16, a steady solution being reached after 5 cycles and used Makie.jl to visualise results [62].

Results
In this section, we present: (i) the convergence of total order indices with respect to both discrete and continuous output measurements; (ii) our investigation outcomes of varying the four estimators with the five sampling methodologies defined in Section 2.2.First, the convergence results will be illustrated in Section 3.1.Then in the next two subsections, total order Sobol indices and orthogonality of input parameters for the 1-chamber,

Convergence and uncertainty
Fig 4 shows the convergence and the uncertainty of the minimal ventricular elastance E min for the 1-chamber and 2-chamber models.We calculate the Sobol indices using the Jansen estimator and Sobol sampling, which are considered to be best practice [18,49].Henceforth, they are regarded as our benchmark and the sample size returned from this initial investigation is used for evaluations on all other estimators and sampling methods.Increasing the sample size and re-sampling with replacement allow us to evaluate the sample size at which the Sobol indices have converged with minimum uncertainty.This is displayed as a band around the index of interest and represents a 95% confidence interval of the index estimate.

1-chamber model
The uncertainty associated with the computation of Sobol total order indices on the 1-chamber model, with N = 10, 000 samples, is presented in Fig 5 .Only the results for arterial compliance C sa are displayed in this figure.The 95% confidence intervals for the Homma and Sobol estimators are considerably wider than that of the Jansen and Janon estimators.The Homma estimator consistently produced estimates of the sensitivity indices which are different to that of the other available estimators.The Jansen and Janon estimators are identical in their computations of the sensitivity indices and confidence intervals, invariant of the sampling methodology used.When the Homma and Sobol estimators are used, the latin hypercube and uniform sampling methods produce larger confidence intervals compared to the quasi-monte carlo sampling methods.
When calculating total order indices on the 1-chamber model with continuous measurements (the histograms for the orthogonality distributions of input parameters are presented in Fig 6), we notice that the orthogonality spreads for the Jansen and Janon estimators are identical for the golden ratio and Latin hypercube sampling methodologies.The Jansen estimator coupled with lattice rule sampling also shares this orthogonality distribution.The Janon estimator, with Lattice Rule and Sobol sampling, and the Jansen estimator with Sobol sampling, exhibit minor variations from the previous orthogonality distributions generated by Janon and Jansen estimators, however, are identical between themselves.The orthogonality distributions returned from the Homma and Sobol estimators exhibit large variations for each sampling methodology, although the Sobol estimator with the Sobol sampling returns an orthogonality distribution similar to that seen by the Jansen and Janon estimators.These results are mirrored in Table 4 where the input parameters are ranked based on their orthogonality scores in the parameter space.We see the orthogonality results obtained for the Jansen and Janon estimators are invariant to sampling methodologies.In contrast, the rankings for Sobol and Homma estimators vary amongst different sampling methodologies.The Sobol estimator when coupled with Sobol sampling, returns a parameter ranking almost identical to that of the Jansen and Janon estimators, a result consistent with the one observed in Fig 6 .In Table 5, stratification by estimator type and examination of the range of an input parameter across all sampling methodologies reveal, as inferred from Table 4, that the Jansen and Janon estimators exhibit no variation for the whole input parameter set, given any sampling methodology.This indicates that the Jansen and Janon estimators are the optimal choices for this model.The Homma and Sobol estimators exhibit variations of 1.33 and 1.67 upon the input parameter set, respectively.These variations mean that using Homma and Sobol estimators will return differing orthogonality rankings when different sampling methodologies are used.When stratifying by sampling types, Table 6 reveals that Sobol and Lattice Rule samplings exhibit the smallest mean variations of the input set across all estimator types.It is important to note that these variations are a consequence of the Sobol and Homma estimators which both exhibited different orthogonality rankings for input parameters.These results indicate that given a less than optimal estimator, the Sobol or Lattice rule sampling methodology may produce a ranking which can be considered closer to the ground "truth".Interestingly, we notice that the commonly used Latin Hypercube sampling methodology in life sciences exhibits the largest variation of an input set of parameters.measurements.In all cases, as the sample sizes are increased, the accuracy of the estimations and uncertainty associated with the indices improve.The Jansen and Janon estimators provide the most efficient convergences and smallest errors when calculating the indices.A sample size of N = 10, 000 is taken for the discrete measurements, because the columns for the Jansen and Janon estimators (Panels K-T) illustrate that any additional samples would return minimal improvements in terms of accurate calculation of the indices.The Homma and Sobol estimators (Panels A-J) display considerably larger errors than that of the Jansen and Janon estimators.When the upper limit sample size of k = 40, 000 is reached, the Homma and Sobol estimators appear to have converged with reduced errors when combined with the Sobol, Lattice Rule and Golden sampling methods, although the errors are still much larger than those exhibited by the Jansen and Janon estimators.The uniform and Latin hypercube sampling methods present the largest errors when combined with the Sobol estimator.
When calculating total order indices on the 1-chamber model with discrete measurements, the histograms presented in Fig 8 show that the orthogonality spreads for the Jansen and Janon estimators for all sampling methodologies, except the Janon estimator and the Latin Hypercube sampling pairing, are identical.We notice, the orthogonality distributions returned from the Homma estimator exhibit large variations for each sampling methodology, as seen with continuous measurements shown in Fig 6 .The Sobol estimator column (Panels F-J) in  7 where the Jansen and Janon estimators are shown to invariant to sampling methodologies, whilst the rankings of parameter orthogonality for the Sobol and Homma estimators vary amongst different sampling methodologies.In Table 8, stratifying by estimator type and examining the range an input parameter exhibits across all sampling methodologies reveals that the Jansen and Janon estimators exhibit no variation for any input parameter set, given any sampling methodology, once again implying they are the optimum choice.The Homma and Sobol estimators exhibit variations of 5.11 to 2.22, respectively, upon the input parameter set.The variations for the discrete measurements are much greater than the variations seen with continuous measurements.When stratifying by sampling type, Table 9 shows the Sobol sampling method exhibits the smallest mean variation of an input set across all estimator types.As above, due to only the Sobol and Homma estimator exhibiting largely varying parameter rankings, stratifying by sampling methodology places more emphasis on the apparent less robust estimators.From this result, it does appear the Sobol sampling method may improve the robustness associated with a parameter orthogonality ranking.One notable observation as seen in Tables 4 and 7 is that while the robustness of the Jansen and Janon estimator can be observed in both tables, the ranking associated with the orthogonality of input parameters changes quite dramatically (for example, E min , R mv , Z ao and R s ), highlighting how the change in data type may have consequences in parameter interpretation when conducting a sensitivity analysis.
Overall, for the 1-chamber model with 9 input parameters, we have seen consistent themes of the Jansen and Janon estimators being most robust and most reliable which appear attributable to the excellent convergence exhibited by these estimators.Sobol and Homma estimators exhibit very variable parameter rankings across different sampling methodologies which are in line with the poor convergence of these estimators.The Sobol sampling method appears to reduce the level of uncertainty associated with an input parameter's orthogonality ranking, as shown in Tables 6 and 9. Interestingly, continuous measurements appear to reduce the level of variation associated with parameter orthogonality ranking when compared to discrete measurements.

2-chamber model
In this section, we present the results of our uncertainty study on the 20 parameters, 2-chamber model.First, the uncertainty associated with the computation of Sobol total order indices, for this more complex model, with N = 20, 000 samples, is presented in Fig 9 .Only the results for the left ventricular maximum elastance E max lv are displayed here.The Jansen and Janon   the Jansen and Janon estimators are nearly identical for all sampling methodologies, excluding the uniform sampling, which shares the same orthogonality spread between the two estimators.With the Homma estimator, the spreads of orthogonality appear more consistent amongst the sampling techniques associated with itself, comparing against the test cases on the 1-chamber model, however, they are largely different from what returned from the Jansen and Janon estimators.The Sobol estimator generated results which appear closer to the orthogonality distributions of the Jansen and Janon estimators, whilst they are not identical, the orthogonality distributions between different sampling methodologies are more consistent than the ones exhibited by the Homma estimator.Examining Table 10, the rankings of input parameters are more consistent for the Jansen and Janon estimators, although there are slight discrepancies, as seen on the simple 1-chamber model.
In Table 11, stratifying by estimator type and examining the range an input parameter exhibits across all sampling methodologies reveal that the Jansen and Janon estimators exhibit minimal variations to sampling methodologies-1.3 and 1.45, respectively.The Homma and Sobol estimators exhibit variations of 8.2 and 7.35 respectively upon the input parameter set.When stratifying by sampling type, Table 12 shows the Sobol sampling method exhibits the smallest mean variation of an input set across all estimator types of 7.2.This is still a large variation due to all estimator types being considered and therefore the range accounts for some of the spurious values generated by the Homma and Sobol estimators.
Fig 11 displays the convergence and uncertainty associated with the computation of the total order indices of the venous compliance C svn , for the 2-chamber model, against the discrete measurements.We see in all cases that as the sample size is increased, the estimate and uncertainty associated with the indices improve.Similar as the continuous measurement case for E max lv shown in Fig 9, the Jansen and Janon estimators provide the most efficient convergence and the smallest error when calculating the indices.A sample size of N = 20, 000 is taken for the discrete measurements, as seen in the Jansen and Janon columns (Panels K-T), any additional sampling would return minimal improvements in terms of accurate calculation of the indices.The Homma and Sobol estimators display errors which are considerably larger than that of Jansen and Janon estimators.We see when the upper limit sample size of N = 30, 000 is reached, the Homma and Sobol estimator errors are still large.This results demonstrates that for this complex model, less efficient estimators (such as Homma and Sobol) and a less accurate sampling method (such as Latin hypercube) display large confidence intervals and struggle to return converged index values.
From the histograms presented in Fig 12, the orthogonality spreads exhibit similar trends to that of the continuous measurements, shown in Fig 10 .We note that the histograms are identical for the Jansen and Janon except when the Uniform sampling method is used (which exhibits slight variations from the other histograms).With the Homma and Sobol estimators, although there appears to be low level consistency amongst their orthogonality distributions, they are very different to the ones produced by the Jansen and Janon estimators.Examining Table 10.Input parameter ranking for the 2-chamber model with continuous measurements-Here, input parameters are ranked based on the averaged orthogonality score returned from the calculated total order sensitivity matrix.In addition, the ranking is stratified by both sampling and estimator types.Table 13, the rankings of input parameters are consistent for the Jansen and Janon estimators apart from the Uniform column which often returns a parameter ranking differing to the other sampling types.

Homma
In Table 14, stratifying by estimator type and examining the range a input parameter exhibits across all sampling methodologies reveals that the Jansen and Janon estimators exhibit minimal variation to sampling methodologies-0.9 and 1.0, respectively.This is an improvement on the continuous measurements as the Jansen estimator returns less than 1 parameter range variation.The Homma and Sobol estimators produce variations of 10.9 and 5.55 respectively upon the input parameter set.When stratifying by sampling type, Table 15 shows the Lattice Rule sampling method exhibits the smallest mean variation of an input set across all estimator types of 5.75.This is still a large variation but is less than the variation exhibited by the continuous measurements.
Overall, for the more complex 2-chamber model with 20 input parameters, the Jansen and Janon estimators are consistently the most robust and reliable estimators.When using continuous measurements, neither returns an input parameter set mean variation greater than 1.When using discrete measurements, they return a mean variation less than or equal to 1 (see Tables 11 and 14).This, as in the 1-chamber case, could be attributable to the efficient rate of convergence displayed by the Jansen and Janon estimator.The Sobol and Homma estimators exhibit very different parameter rankings across different sampling methodologies with variations of up to 10.9.These large variations are in line with the poor convergence associated with these estimators.The Sobol and Lattice Rule sampling method appears to reduce the level of uncertainty associated with an input parameter's orthogonality ranking (see Tables 12 and 15), despite spurious parameter rankings from the Homma and Sobol estimators leading to large parameter variation when stratified by sampling methodologies.

Discussion
Utilising two cardiovascular system models, the main aim of our investigation is to test the robustness of the calculation of the input parameter orthogonality, while varying total order estimator types and sampling methodologies, across differing input parameter dimensionalities and types of data on which the total order indices were calculated.The results presented in Section 3 display overwhelming robustness for the Jansen and Janon estimators when calculating total order indices when compared to other options.For the 1-chamber, 9-parameter model, we observed that these two estimators gave nearly invariant outcomes to both the sampling methodology and the data type.When the dimensionality of the model parameters is increased to 20, we noted that the Jansen and Janon estimators exhibited small variations on the input parameter orthogonality rankings.For the Jansen estimator with discrete measurements, it returned a mean variation of less than 1.We observed that the Homma and Sobol estimators regularly returned mean variations for input parameter sets greater than 1, which is particularly amplified when the model dimensionality is increased.
Given our aim is to assess the use-ability of estimators and sampling methodologies for practical identifiability studies, these results indicate that if the used estimator and sampling methodology are not robust, the calculated optimal parameter set is unreliable.The usage of unrobust estimators and sampling methods can therefore produce misleading conclusions with practical consequences, especially in investigations applied to life sciences.Interestingly, we witnessed that the variations attached to the Sobol and Lattice Rule sampling methods were the lowest across all model dimensionalities and data types.Our results also reinforce the findings reported in [14,63] that the commonly used Latin Hypercube method is less than optimal in exploring the input parameter space, especially at high dimensionalites.The Jansen and Janon robustness at calculating total order indices can be attributed to the Jansen estimator never allowing negative values in the numerator (see estimator definitions in Table 1), where as the Janon estimator is the only estimator which has been proven to be both asymptotically normally distributed and asymptotically efficient meaning as the sample size increases the estimation error associated with calculating the indices is negligible [50].They are both highly optimised estimators with very little room for improvement [49,50].The allowance of negative indices in both the Sobol and Homma estimators is an explanation for the poor performance in the calculations of these indices.From the Hoeffding-Sobol decomposition, the total order indices can be viewed as a decomposition of the variance.Thus, any negative value obtained from an estimator is not theoretically possible (the average of the squared deviations).Any negative value returned from an estimator is likely to be due to numerically instability, or some points in the parameter space returning unphysiological results.Thus reinforcing the important point of ensuring any sensitivity indices have adequately converged in order to avoid any of these unstable points.
While the robustness of the Homma and Sobol estimators have been questioned within this work, this does not mean they lack ability to produce accurate estimations of parameter influence.Fig 13 displays the Homma and Sobol total order estimator results against the maximum left ventricular volume with an increased sample size comparing to the tests conducted earlier -N = 100, 000 (which requires 2, 200, 000 model evaluations).We observe that now with a much extended sample size, the confidence intervals (green bands in Fig 13) and parameter   interpretations from Homma and Sobol estimators are similar to that of the Jansen or Janon estimators, which are displayed in panels C and D for N = 40, 000 samples.Theoretically, this demonstrates all estimators can produce equivalent results, when 'large enough' samples are used.However, in the case of a limited computational budget, our results indicate a clear estimator choice, to obtain robust and valid parameter interpretations with the best computational efficiency.
With the highly optimised Jansen and Janon estimators, it is evident that they consistently exhibit the most efficient convergence and produce the smallest uncertainties when calculating total order indices.This phenomenon matches the consistent orthogonality observed among input parameters across various sampling techniques when coupled with the Jansen and Janon estimator.Conversely, the Homma and Sobol estimators tend to yield significantly larger uncertainties when sample sizes are held constant among estimators, thus explaining the lack of consistent orthogonality rankings for the input parameters.Increasing the sample size seems to ameliorate the uncertainties associated with the Homma and Sobol estimators, particularly when employing the Sobol, lattice rule, and Golden sampling methods.This observation underscores the resilience of low-discrepancy sequences, demonstrating their effectiveness even in conjunction with a sub-optimal estimator.While the work of Puy et al. [18] did not delve into extensive convergence or uncertainty quantification, it is plausible to infer that the Jansen and Janon estimators, with their superior convergence rates and lower uncertainty, played a pivotal role in their conclusion that these estimators are the most efficient at capturing the true effects of input parameters.
Our investigation has been confined to two highly non-linear stiff differential algebraic equation systems with the understanding that they represent a high level of complexity, therefore good modelling guidelines obtained here would be readily applicable to simpler, more linear models which are associated with a less variable input parameter space.Thus, obtaining sensitivity estimates for linear models is considerably less expensive than what has been conducted here.Convergence and uncertainty quantification have historically been left out of sensitivity analysis studies, despite being highlighted as vital, if the results of studies were then to influence policy /societal /clinical decisions [14,64].In this study, we have highlighted the impact that convergence has on a total order estimator, alongside this, we have also shown the impact of the level of sampling taken, on the calculation of total order estimators.It appears intuitively sensible that the higher density of sampling leads to better resolution of the input parameter space, hence our sensitivity analysis gives a better indication about which input parameters are truly influential.Current literature states that N > 500 but this recommendation is based on physical systems which are mostly linear.The work conducted and results shown in this study highlight that for a highly non-linear system, one should investigate N > 5000.More importantly, it is clear that no two systems are the same, so for one to ensure adequate resolution of an input parameter space, convergence and uncertainty quantification through bootstrapping must be an essential step in any modeller's workflow, given the aim is to perform accurate and robust parameter identification studies.
While it is clear that a large sample size is needed in order to ensure accurate sensitivity estimates, researchers are often restricted by the computational budget available to them.Without sufficient computational speed, it would be unfeasible to perform the needed number of model evaluations in order to obtain accurate indices interpretation.In our example of the 2-chamber model, 100k samples required 2.2 million model evaluations.1 model evaluation took 0.039 seconds to run in Julia, thus would require 23.8 hours if this was computed in serial.The Julia language has consistently shown to be 10 to 100 times [59,65] quicker than other comparable languages in solving systems of this class and the usage of parallel computation on high performance computing facilities has further sped up our simulations.Therefore, in situations of a limited computational budget, the selection of an efficient language and estimators which require less samples to reach convergence is key.
Although not exhaustive, the list of estimators investigated in this paper does represent what are readily available, practically usable and computationally feasible.Puy et al. [18] also recommended the estimator introduced by Azzini et al. [66], which appeared to produce similar results to the Jansen estimator.The Azzini estimator requires k = 2N(n + 1) model evaluations, compared to N(n + 2) evaluations needed by the estimators investigated in this work.The 2-chamber model would require k = 1, 320, 000 model evaluations for Azzini estimator.For models with higher parameter dimensions, this would be computationally infeasible.With the increasing prospect of digital twins in healthcare, more complex and detailed models which accurately represent the true physiological processes are generated.However, the bottleneck which prevents the progression of these models to clinically applicable situations is the computational cost associated with a detailed sensitivity analysis.This does not refer to the calculation of the indices, but the process of solving the dynamical system.Therefore, while new estimators may prove to be accurate, the focus must be on efficient resolution of complex dynamical systems and the efficiency of the estimators for low sample numbers, in order to ensure a thorough sensitivity analysis.
All the estimators used in this work are available in the global sensitivity packages such as SALib, GlobalSensitivity.jl,SenSobol and sbiosobol [60,67,68].It is reassuring to see that the default estimator used to calculate the total order index, in the available packages, is the Jansen estimator.Given the conclusion drawn from Puy et al. [18] and the findings from this work, researchers could straightforwardly use the above packages when performing practical identifiability studies and would obtain a reliable optimal set of input parameters which best describes the experimental data available to them.
Another area of research is the calculation of total order sensitivity indices where one assumes dependency between input parameters.This is partially investigated by Puy et al. [68] in implementing the method of Glen et al. [69].This method requires a prescription of linear dependencies between parameters, however, these are often not known in realistic life science models.As a result, Puy et.al. [68] demonstrated the inaccuracy associated with this method when calculating the true effects.There have been various methods deriving variance based sensitivity indices with dependent inputs [70][71][72], however, similar as the method of Glen et al. [69], they require knowledge of the dependencies that exist within the model and therefore the computational power needed to simulate these indices is often much larger than the standard Sobol indices.On top of this, there is no accepted method for how to calculate these dependent indices which should be of interest for future work.Therefore, the need to understand how input parameter orthogonality is affected by varying estimators and sampling methodologies is of significant importance, in order for total order sensitivity indices to be utilised in identifiability studies.
While we have conducted this work through the lens of utilising the method of Li et.al. [7] (see Eq (1)) for practical identifiability studies, it is also applicable to other methods.Another approach of identifying input parameters is the structured correlations method [73], where one seeks to identify correlations between parameters and to calculate ranks (based on which parameters can be identified uniquely if they are not strongly correlated with other parameters).This approach utilises the total order sensitivity matrix to calculate these correlations.Therefore, the need for reliable and robust sensitivity matrices is vital to whichever method is implemented.
The work of Puy et al. [18] is conclusive in its findings of the Jansen and Janon estimators being the most reliable in finding the "truth" input parameter effects.Our work complements these conclusions in that we find the Jansen and Janon are the most reliable in the calculation of input parameter orthogonality, which appears to be motivated by input parameter convergence.Puy et al. [18] also reported that any choice made on the model have a non-negligible effect.While we agree with these conclusions for the most part, we are able to identify that similar to the model dimensionality, it appears choices we make, such as the sampling methodology and type of data used, have more impactful consequences.The lower variation associated with input parameter sets, when low discrepancy sequences are used, implies their effectiveness in returning robust and reliable input parameter sets.It appears that there is no clear advantage of using either continuous or discrete measurements when choosing how to calculate the total order indices.

Conclusion
Our study delved into the intricacies of varying sampling methodologies and variance-based total order estimators, aiming to establish best practices for practical identifiability studies.We conducted our investigation using two highly non-linear and stiff 0D models of the human cardiovascular system as our test cases: (i) a 1-chamber, 9-parameter model and (ii) a 2-chamber, 20-parameter model, both based on differential algebraic equations.Through a thorough empirical assessment of total order estimators and sampling methodologies, we gained valuable insights into their strengths and weaknesses, shedding light on the orthogonality of the input parameters within the models.This analysis complements prior work that focused on the estimators' ability to uncover the "true" effects of a model, enriching our comprehension of their practical identification.
Our findings strongly advocate for the Jansen and Janon estimators as robust choices across different sampling methodologies, measurement data variations, and model dimensions.These two estimators emerge as preferred tools for calculating total order indices and, consequently, for identifying the optimal set of input parameters.Their efficient convergence and the consequential reduction in index uncertainty make them the optimal choice for this task.Furthermore, we recommend the use of low-discrepancy quasi-random Sobol and Lattice Rule sampling schemes as optimal sampling methodologies to complement Jansen and Janon estimators.
In essence, our work establishes a robust framework of good modelling practice for practical identifiability studies, considering both the influence of input parameters and their orthogonality.By incorporating these best practices into modeling studies, researchers can consistently and reliably identify the optimal input parameters for dynamical systems.This approach not only enhances the quality and accuracy of parameter identification, but also paves the way for more informed decision-making in various scientific and practical domains.
Fig 2 details the workflow needed to generate the sensitivity indices in Julia.
9-parameter model (Fig 1A) and the 2-chamber, 20-parameter model (Fig 1B) are shown.Between each subsection, we examine what effect the different choices of estimator and sampling methodology has on the orthogonality of input parameters.Within each of the subsections, we make the distinction between the effects of continuous and discrete measurements.We present results for input parameters which are deemed to have high clinical significance (i.e., biomarkers), for example, low arterial compliance C sa may indicate a stiffening of the vessel.Each subsection displays convergence of a single parameter for brevity.All other data are available at https://github.com/H-Sax/Orthgonality-SA.Physiologically realistic time series solutions generated from the model is shown in Fig 3.
Fig 4  shows the convergence and the uncertainty of the minimal ventricular elastance E min for the 1-chamber and 2-chamber models.We calculate the Sobol indices using the Jansen estimator and Sobol sampling, which are considered to be best practice[18,49].Henceforth, they are regarded as our benchmark and the sample size returned from this initial investigation is used for evaluations on all other estimators and sampling methods.Increasing the sample size and re-sampling with replacement allow us to evaluate the sample size at which the Sobol indices have converged with minimum uncertainty.This is displayed as a band around the index of interest and represents a 95% confidence interval of the index estimate.For the 1-chamber model, 10, 000 samples (110, 000 model evaluations) ensured convergence when computed against the discrete measurements defined in Eq(28).Fig4Ashows that evaluating the Sobol indices at a higher sample size would provide minimal improvement

Fig 3 .
Fig 3. Time series solutions for the 1 and 2 chamber cardiovascular models investigated in this work.The solutions shown are the ones which are utilised in the investigation.https://doi.org/10.1371/journal.pcbi.1011946.g003

Fig 4 .
Fig 4. Convergence and uncertainty of indices associated with the minimum ventricular elastance E min .Fig A displays the convergence and uncertainty of the Sobol indices S T calculated on discrete measurements for the 1-chamber model against increasing sample size.Here, the vertical line signifies the chosen sample size for the 1-chamber model at N = 10, 000.Fig B presents the continuous Sobol indices with uncertainty bounds, calculated at a sample size N = 10, 000, on continuous measurements over a single cardiac cycle, for the 1-chamber model.Fig C displays the convergence and uncertainty of S T calculated on discrete measurements for the 2-chamber model against increasing sample size.Again, the vertical line signifies the chosen sample size for this model, at N = 20, 000.Fig D shows the continuous Sobol indices with uncertainty bounds for N = 20, 000, on continuous measurements over a single cardiac cycle, for the 2-chamber model.The measurements shown in blue, yellow and green denote the left ventricular pressure, the systemic arterial pressure and the left ventricular volume, respectively.In the discrete settings (i.e., A and C), the measurements are the mean left ventricular pressure, the maximum systemic arterial pressure and the maximum left ventricular volume.https://doi.org/10.1371/journal.pcbi.1011946.g004 Fig 7 displays the convergence and uncertainty associated with the computation of the total order indices of the mitral value resistance R mv for the 1-chamber model against the discrete

Fig 5 .
Fig 5. Total order Sobol indices S T of the arterial compliance C sa for the 1-chamber model with continuous measurements.Panels A-T show S T of C sa , for 3 continuous measurements-left ventricular pressure, systemic arterial pressure and the left ventricular volume (represented in blue, yellow and green curves, respectively), over a single cardiac cycle with differing estimators and sampling methodologies.Measurements are evaluated with N = 10, 000 samples, using B = 1000 bootstrapped samples to evaluate the uncertainty of the estimate.The bands represent 95% confidence intervals associated with specific indices displayed as solid curves.https://doi.org/10.1371/journal.pcbi.1011946.g005

Fig 6 .
Fig 6.Orthogonality distributions of input parameters for the 1-chamber model with continuous measurements-Histograms A-T show the distribution of orthogonality returned from examinations of the sensitivity vectors, calculated from continuous measurements.Here, an orthogonality score of 1 represents total independence of input parameters, whereas 0 represents total dependence.Each individual diagram denotes a specific combination of sampling methodology and estimator type.The frequency of each histogram is normalised such that it is comparable between plots, i.e., the larger the frequency of a bin, the larger the number of orthogonality scores calculated from the original sensitivity vectors.https://doi.org/10.1371/journal.pcbi.1011946.g006

Fig 7 .
Fig 7. Total order Sobol indices S T of the mitral valve resistance R mv for the 1-chamber model with discrete measurements.Panels A-T show S T of R mv , for 3 discrete measurements: mean left ventricular pressure, maximum systemic arterial pressure and maximum left ventricular volume (represented in blue, yellow and green, respectively), evaluated at increasing sample sizes (N 2 [2000, 40000] using B = 1000 bootstrapped samples), with differing estimators and sampling methodologies.The bands represent 95% confidence intervals associated with specific indices displayed as solid curves.The red solid vertical lines represent the point (N = 10, 000) at which the sample size is taken.https://doi.org/10.1371/journal.pcbi.1011946.g007

Fig 8 .
Fig 8. Orthogonality distributions of input parameters for the 1-chamber model with discrete measurements-Histograms A-T show the distribution of orthogonality returned from examinations of the sensitivityvectors, calculated from continuous measurements.Here, an orthogonality score of 1 represents total independence of input parameters, whereas 0 represents total dependence.Each individual diagram denotes a specific combination of sampling methodology and estimator type.The frequency of each histogram is normalised such that it is comparable between plots, i.e., the larger the frequency of a bin, the larger the number of orthogonality scores calculated from the original sensitivity vectors.https://doi.org/10.1371/journal.pcbi.1011946.g008 estimators display identical index and confidence interval calculations, apart from when combined with the lattice rule sampling methodology, the plots show slightly larger confidence intervals for the left ventricular volume comparing to the other cases.The Homma and Sobol estimators again display much larger confidence interval estimates compared to the Jansen and Janon estimators.The errors associated with the Sobol sampling methodology when the Homma and Sobol estimator are used, are much smaller compared to the Latin hypercube and uniform sampling methodologies, hence demonstrating the impact sampling methodology can have on the estimations of sensitivity indices.Next, we calculate total order indices on the 20 dimensional 2-chamber model with continuous measurements.From the histograms presented in Fig 10, the orthogonality spreads for

Fig 9 .
Fig 9. Total order Sobol indices S T of the maximal left ventricular elastance E max lv for the 2-chamber model with continuous measurements.Panels A-T show S T of E max lv , for 3 continuous measurements-left ventricular pressure, systemic arterial pressure and the left ventricular volume (represented in blue, yellow and green curves, respectively), over a single cardiac cycle with differing estimators and sampling methodologies.Measurements are evaluated with N = 20, 000 samples, using B = 1000 bootstrapped samples to evaluate the uncertainty of the estimate.The bands represent 95% confidence intervals associated with specific indices displayed as solid curves.https://doi.org/10.1371/journal.pcbi.1011946.g009

Fig 10 .
Fig 10.Orthogonality distributions of input parameters for the 2-chamber model with continuous measurements-Histograms A-T show the distribution of orthogonality returned from examinations of the sensitivity vectors, calculated from continuous measurements.Here, an orthogonality score of 1 represents total independence of input parameters, whereas 0 represents total dependence.Each individual diagram denotes a specific combination of sampling methodology and estimator type.The frequency of each histogram is normalised such that it is comparable between plots, i.e., the larger the frequency of a bin, the larger the number of orthogonality scores calculated from the original sensitivity vectors.https://doi.org/10.1371/journal.pcbi.1011946.g010

Fig 11 .
Fig 11.Total order Sobol indices S T of the venous compliance C svn for the 2-chamber model with discrete measurements.Panels A-T show S T of C svn , for 3 discrete measurements: mean left ventricular pressure, maximum systemic arterial pressure and maximum left ventricular volume (represented in blue, yellow and green, respectively), evaluated at increasing sample sizes (N 2 [10000, 30000] using B = 1000 bootstrapped samples), with differing estimators and sampling methodologies.The bands represent 95% confidence intervals associated with specific indices displayed as solid curves.The red solid vertical lines represent the point (N = 20, 000) at which the sample size is taken.https://doi.org/10.1371/journal.pcbi.1011946.g011

Fig 12 .
Fig 12. Orthogonality distributions of input parameters for the 2-chamber model with discrete measurements-Histograms A-T show the distribution of orthogonality returned from examinations of the sensitivity vectors, calculated from continuous measurements.Here, an orthogonality score of 1 represents total independence of input parameters, whereas 0 represents total dependence.Each individual diagram denotes a specific combination of sampling methodology and estimator type.The frequency of each histogram is normalised such that it is comparable between plots, i.e., the larger the frequency of a bin, the larger the number of orthogonality scores calculated from the original sensitivity vectors.https://doi.org/10.1371/journal.pcbi.1011946.g012

Fig 13 .
Fig 13.Estimator comparisons with larger samples for the 2-chamber model with discrete measurements-The Sobol and Homma estimators results are based on 100k samples, compared to the Jansen and Janon estimators using 40k samples both with 95% confidence.The input parameter effect is displayed against the maximum left ventricular volume as an example here.https://doi.org/10.1371/journal.pcbi.1011946.g013

Table 2 . Input parameters for the 1 chamber model.
Each input parameter's unit is stated alongside a chosen initial value for the 9 parameter, 1-chamber model.τ is the cardiac cycle length and is fixed such that τ = 1s.The ventricular shift parameter E shift = 0 s as no atrium is present in this model.

Table 3 . Input parameters for the 2 chambers model
https://doi.org/10.1371/journal.pcbi.1011946.t003both models are passive, having fixed compliances; (iv) flow in and out of the active left atrium in Fig 1B is controlled by the systemic veins and the mitral valve.For both models, in Fig 1A and 1B, flow in and out of the active left ventricle is controlled by the mitral and aortic valves respectively.The valves are modelled as diodes, with Ohmic resistance under forward bias and infinite resistance under reverse bias: . Each input parameter's unit is stated alongside a chosen initial value for the 20 parameter, 2-chamber model.τ is the cardiac cycle length and is fixed such that τ = 1s.The ventricular shift parameter E shift = 0.92 s as an atrium is present in this advanced 20 parameters model.

Table 4 . Input parameter ranking for the 1-chamber model with continuous measurements-Here, input parameters are ranked based on the averaged orthogonal- ity score returned from the calculated total order sensitivity matrix.
In addition, the ranking is stratified by both sampling and estimator types. https://doi.org/10.1371/journal.pcbi.1011946.t004

Table 7 . Input parameter ranking for the 1-chamber model with discrete measurements-Again, input parameters are ranked based on the averaged orthogonality score returned from the calculated total order sensitivity matrix.
The ranking is also stratified by both sampling and estimator types. https://doi.org/10.1371/journal.pcbi.1011946.t007

Table 12 . The ranges of input parameters across 4 estimator types for a specific sampling method for the 2-chamber model with continuous measurements.
https://doi.org/10.1371/journal.pcbi.1011946.t012

Table 13 . Input parameter ranking for the 2-chamber model with discrete measurements-Here, input parameters are ranked based on the averaged orthogonality score returned from the calculated total order sensitivity matrix.
In addition, the ranking is stratified by both sampling and estimator types.

Table 14 . The ranges of input parameters across 5 sampling types for a specific estimator for the systemic circulation model with discrete measurements.
https://doi.org/10.1371/journal.pcbi.1011946.t014

Table 15 . The ranges of input parameters across 4 estimator types for a specific sampling method for the 2-chamber model with discrete measurements.
https://doi.org/10.1371/journal.pcbi.1011946.t015