Assessing input parameter hyperspace and parameter identifiability in a cardiovascular system model via sensitivity analysis

We aim to clarify our understanding of the process of state-space model input parameter identification, known, within the clinical context, as model personalisation. To do so, we apply reference sensitivity and identifiability techniques to a lumped parameter, single ventricle representation of the systemic circulation, chosen in view of its relative simplicity and prior art. We attempt to quantify the reliability of input parameter identifiability through the lens of 4 clinically relevant measurements and the attendant difficulty in personalising the model. In turn, this that we extend existing methods which combine both parameter influence and orthogonality, to global sensitivities. By examining different parameter sensitivity evaluation methodologies, we investigate the stability of optimal parameter subsets which are commonly used to aid clinical investigations. In order to perform the personalisation process, one must understand the complexity of the high dimensional input parameter hyperspace associated with this class of model. By utilising Sobol indices, we propose a domain-agnostic and intuitive approach. This involves varying the bounds of the input parameter space relative to the model’s base state. These investigations yield a pseudo-mapping of the input hyperspace, cementing our understanding of the role of identifiable input parameters in the state-space model. Our findings suggest a novel global methodology for input parameter identifiability and input hyperspace mapping, providing valuable insights into solving the personalisation process.


Introduction
Dynamical systems in the life sciences are mathematical descriptions of processes observed in nature [1].These systems can be used to predict and infer a wide range of biological processes, such as population dynamics [2], cell mechanics [3] and global hemodynamics [4].Due to the complexity of biological processes, the corresponding dynamical systems are often described through a large set of nonlinear equations, which take the form of ordinary differential equations (ODEs), differential algebraic equations (DAEs), partial differential equations (PDEs) or a mixture [5,6].Due to complexity of the formulation, computational tools are required to solve the model and make inferences about the processes under investigation.An invariant property between all dynamical systems is that they are parameterised by the set of values of input parameters, which themselves can provide insight into the process being modelled [7] (e.g., species carrying capacity, drug diffusion rate and systemic vascular resistance).
One impactful use of dynamical systems, and the focus of this study, is the modelling of the cardiovascular system (CVS).CVS models are often categorised by their dimensionality, each class of models having a very different purpose.Zero-dimensional (0D), or lumped parameter models, divide the CVS into compartments within which the state variables are assumed to be uniformly distributed and vary only with time.0D models can be used to represent the whole CVS physiology, or any portion of it [8].The physiological state variables of pressure, flow and volume are respectively equivalent to voltage, current and charge in an electrical analogue.Each CVS model or compartment can be represented as a combination of resistors, capacitors and inductances which are parameterised by numerical values ,  and  respectively.For a generic vessel or organ located in a larger circulation network, ,  and  represent hemodynamic dissipation, vessel distensibility and the inertial effects blood flow, respectively [9].Thus, the input parameters of 0D models carry great clinical significance as bio-markers, used https://doi.org/10.1016/j.jocs.2024.102287Received 30 October 2023; Received in revised form 28 March 2024; Accepted 3 April 2024 to inform medical treatment [10].Higher dimensional models of the CVS which include spatial variation, are also investigated.This class of model is based with pdes and is usually utilised for detailed investigation of specific vessels within the wider circulation network [11][12][13][14], owing to the higher computational cost associated with solving PDEs, rather than ODEs/DAEs.In order to tackle the computational cost reduced order modelling is a suitable choice [15,16].Also, the input parameter list is often larger and contains less clinically relevant parameters, making output data model analysis harder.Therefore, 0D models of the CVS are our modelling technique of choice in this work.For a detailed review of modelling types see e.g.[17,18].
Mathematically, input parameter identification is an inverse problem; clinically it is the personalisation process [19,20].It amounts to taking a highly detailed dynamical system model with its input parameters (which are prescribed realistic bounds) and supplying the model with target observations i.e. experimental or clinical data that correspond to model outputs.Performing the process returns a new set of input parameters within the known, realistic bounds, which best describe the observed data.This new set allows for new inferences to be made about the physical process under investigation [21] or the condition of the patient from whom the target data was taken.In the clinic, 0D models have shown great promise in aiding the personalised diagnosis and treatment of CV diseases such as coronary artery disease, pulmonary arterial hypertension and aortic valve stenosis [22][23][24].If one can inform CVS 0D models with measurements specific to patients and personalise quickly and robustly, one may hope to remove the need for invasive diagnostic tests.The diagnostic challenge, then, requires that a model be able to personalise to patho-physiological states, as well the physiological.Further, if personalised input parameters can be used to e.g.stratify a cohort, 0D models might rise still further, to meet the prognostic challenge [25].
The mathematical essence of the personalisation process is a solution for an input parameter coordinate set that locates the global minimum of a hypersurface, or landscape, spanned by the input parameters and computed from the target clinical measurements.Patientspecific data are sparsely available and mathematically insufficient, so many off-line investigations must be performed to ensure the optimal accuracy, efficiency and uniqueness of the solution to the personalisation process.Despite progress there remain many open questions surrounding the key personalisation issues, which we distill as questions • What is the most effective and stable methods for mapping the bounded, physiologically realistic input parameter space?• How does one ensure biomarkers extracted from input parameters are truly patient specific?• What is the surface complexity of input parameter space corresponding to the available measurements?
This study aims to provide methodologies and limited answers for all three questions above.It is important to note our investigation uses forward data, generated from the prescribed dynamical system, in order to understand the methods in an ideal setting.Of course, without this critical, off-line step, Misleading results may be obtained which would then lead to ill-informed clinical decisions.
In Section 2, we review relevant literature, introduce concepts germane to the personalisation process, detail the position this type of investigation has in the personalisation process and summarise the principal contributions of this work.In Section 3, the mathematical detail is provided, for each tool used for model analysis (Sobol indices, identifiability analysis and hyperspace investigation algorithm) and we detail the quantitative investigation of the complexity of input parameter space.Section 4 declares our results from different computational experiments.Discussion of the personalisation process and of the results is given in Section 5.
The key problem is how to effectively navigate the input parameter hyperspace.From the discussion in Section 1, it prudent to review terminologies, prior art and state how the contributions of this work will fresh perspective on the personalisation process.

Terminologies
For the value of input parameters to be identified uniquely from given data, there prerequisite properties.
First, an input parameter effect needs to be influential on the output hyper-surface, i.e., a change in the input parameter space causes a detectable change on the desired output.If one can identify such an input parameter, it is said to be sensitive [29].One can identify locally and globally sensitive input parameters, with respect to the measurements.Locally sensitive input parameters are those eliciting the greatest rate of change of the output about the model base operating point or state [30].Globally sensitive input parameters are potential biomarkers which operate within a physiologically realistic value range.Input parameters are said to be globally sensitive when they cause the greatest influence on the outputs, for the prescribed input range [31].Different methods exist to calculate the sensitivity of input parameters (see Section 3.2).Through the lens of the personalisation process, sensitive input parameters play a vital role.We can first identify the sensitive input parameters -their effect should be present in dataso under inverse operation, the model should be able to identify the change in the sensitive input parameter and possible stratify a patient, with clear clinical implications.
Second, sensitive input parameters having been identified, input parameter orthogonality is considered.Consider two independent influential input parameters with respect to a specific measurement, if the effects of the two input parameters cause the same, or very similar, change in the outputs, the effects of each parameter cannot be extracted individually from the data, under inverse model operation and it is mathematically impossible to identify which input parameters contribute to the effect on the output.It is important for two independent input parameters to have orthogonal effects on specific outputs.When personalising a 0D model, the property of sensitivity is not sufficient, input parameters are also required to be orthogonal.Input parameters that are both sensitive and orthogonal are prime candidates for personalisation [32].Mathematical descriptions of orthogonality are given in Section 3.4.
Personalisation has come to mean the search for an identifiable model and identifiable input parameters.Identifiability analysis of a CVS model needs to include three analyses: structural, sensitivitybased and practical identifiability.Structural (theoretical) identifiability [33] assumes copious and noise free target output data so a model's structural identifiability is deemed clinically academic often.(This, of course, overlooks the fact that one may be unable to identify input parameters because of the model structure, not because of data issues.Naturally, if a model is not structurally identifiable, any practical attempt at its use is limited.)Sensitivity based identifiability analysis [34] combines the above points of identifying sensitive and orthogonal input parameters, subject to model synthetic data, in order to investigate which input parameters are identifiable, in an ideal setting.Practical identifiability analysis [35] accounts for data quality in patient data, where the noise and sampling rate of the patient data may impact the identification of unique input parameters.In the personalisation process, each stage must be performed in order for a model to have direct clinical utility.This study concerns itself with sensitivity based identifiability analysis.

Relevant literature
Our CVS model, (Fig. 1) chosen in part for its relevant prior art, is essentially that proposed by Bjordalsbakke et al. [36].It is a closedloop, lumped parameter model of the systemic circulation and the left ventricle.It is a reasonable candidate for clinical applications, in which focal patho-physiology (say, a defective heart valve) precipitates effects across the patient physiological envelope.Similar models have been proposed, some of which are open-source [37][38][39][40].Bjordalsbakke et al. use the model in Fig. 1 (with a reduced double Hill elastance parameterisation) to demonstrate their Stepwise Subset Reduction Method, or SSRM, for selecting the optimal number of input parameters for the personalisation, using global sensitivity analysis.Bjordalsbakke was also able to confirm the structural identifiability of the model.Using the same essential model, Schiavazzi et al. [38] used Bayesian analysis to examine the probability of recovering unique model parameters and Pironet et al. [41] examined whether model input parameters could be recovered in principle, from chosen outputs.Significantly, Pironet was able to prove injectivity between the inputs and outputs, confirming the structural identifiability of their model.Furthermore, Pironet et al. [42] examined the chances of recovering unique model parameters, when data were corrupted by noise, or were too sparsely collected.Olsen et al. [32] proposed, then examined their structural correlations and orthogonal sensitivities method in reference to a full circulatory model.Whilst a number of other investigations have been conducted, the above cover the main aspects of the personalisation process.Most studies which are concerned with input parameter subset selection do not first understand the identifiability of input parameters in an ideal situation, either by examining just structural considerations or using noisy patient data, thus any identifiability results are often confounded by the innate noise associated with clinical data.Using measurements generated from a model (identifiability in an ideal setting) presents as a natural first step to understand what identifiability issues can be related to data noise.Other limitations reported studies are that local sensitivities are used to select the identifiable input parameters (despite the intrinsic global nature of the personalisation process) and optimisation of a cost function, to some data, without the guarantee that the resulting parameters are unique.
Often, a lack of model sensitivity is deemed the root cause of input parameters' lack of identifiability [38,43,44].Certainly, parameters may be influential but if a group are linearly-dependent (meaning they have similar effects across outputs), it is impossible to determine each uniquely.Accordingly, if one does not examine and screen culprit parameters, obstinately selecting on influence alone, one does not know if a chosen subset of parameters are individually identifiable.A method that aggregates orthogonality and influence is the orthogonal sensitivities method of Li et al. [45], which provides an intuitive means of combining input parameter identifiability with orthogonality.For a full mathematical description see Section 3.5.The original method of Li et al. uses local sensitivities; we argue that an extended exploration of input parameter space is necessitated by the non-local nature of the personalisation problem.Similarly motivated, Ottesen et al. [46] have devised the structural correlations method, which investigates the correlations between parameters of a chosen dynamical system.

Study justification
The personalisation process is a large and detailed procedure with many choices to be made to ensure unique and patient specific input parameters.Given that cohorts of patients may exist at many different patho-physiological states, one must explore the input parameter space globally [21,25].Off-line assessment of a 0D model is a vital prelude.as it ensures a complete understanding, in an ideal setting, of input parameter space.Thus, when practical identifiability is performed, profile likelihood or global optimisation, or any problem in parameter identification can be attributed to issues surrounding available data.
We aim to extend the sensitivity-based identifiability stage of the procedure, to utilise global sensitivities, and to provide an investigative test to quantify the complexity of input parameter, which could aid future studies on practical identifiability.The essential contributions of the work are as follows: Collectively, these contributions enhance the understanding and robustness of the personalisation process, offering valuable tools for parameter identification in patient-specific models.

Methods
Here, we examine the methods used for analysing a declared, minimal system model.In Section 3.1, we declare the model, chosen synthetic measurements and our computational framework.We detail the methods of local and global sensitivity analysis, along with the Fisher information matrix, including the Sobol [47] and eFAST [48] methodologies, in Section 3.2.The definition of orthogonal input parameters given in Section 3.4 and the definition of input parameter influence in Section 3.3.A full description of the global parameter subset selection method in Section 3.5 and in Section 3.6 we detail an investigative test utilising Sobol indices to understand the complexity of the input parameter hyper-surface.

Model and measurements
Mathematically, our system is conveniently expressed in state-space form in which  denotes an input parameter vector,  represents the set of state variables of the system,  is a function describing the system, normally this is an collection of ODEs, ℎ is the measurement function where forward model synthetic measurements are generated, using the computed state variables  and  , as the measurements of interest.

Single ventricle model
The model declared in electrical analogue form in Fig. 1 is a system-level, ordinary differential equation based, electrical analogue CV model, after Bjordalsbakke et al. [36], with three compartments.The state of each compartment is specified by its time-dependent dynamic pressure  (mmHg), inlet flow  (mL/s) and volume  (mL): where  denotes the left ventricle,  the systemic arteries and  the venous system.Formally,  is a continuous variable.
In generic form, the equations relating to the passive compartmental state variables all take the form: Above, the subscripts ( − 1), , ( + 1) respectively represent the proximal, present and distal system compartments,  , (mL) denotes the circulating (stressed) volume [50] and   (ml/mmHg) and   Fig. 1.The systemic circulation, single ventricle, model.An electrical analogue representation of our state space system CVS model, with the elastance of the left ventricle,   , represented by the Shi double cosine model [49].The valves (diodes) are assumed to have Ohmic behaviour, under both forward and reverse bias, with the regurgitating resistance very large.Our textual notation for the resistances and the capacitance's is defined in Table 1.

]
Venous compliance 11.0 (mmHgs/mL) denote compartmental compliance and the Ohmic resistance between compartments , ( + 1).See Fig. 1 and Table 1.Note in our model formulation, (i) we use a C-R-C Windkessel [51] to represent the systemic circulation, (ii) no inertance appears in our formulation, (iii) the systemic and venous compartments are passive, having fixed compliances   and   , respectively and (iv) flow in and out of the active left ventricle is controlled by the mitral and aortic valves, respectively, the latter being modelled as diodes, with Ohmic resistance under forward bias and infinite resistance under reverse bias: where   and   represent the resistance across the aortic and mitral valves respectively.Let us consider the single active model compartment.The dynamics of the left ventricle is described by a time-varying compliance   (), or reciprocal elastance,   () (mmHg/ml), which determines the change in pressure for a given change in the volume [50]: where  0 &  () represent the chamber unstressed and stressed volumes, respectively.() is written [49]: where (;   ,   ) is the activation function for the ventricle and is parameterised by the end systolic and end pulse timing parameters   and   respectively.The elastance function is defined over one cardiac cycle, i.e., time t ∈ [0, ] with  (the length of the cardiac cycle) fixed in this work to Table 2 Derived, discrete output metrics.We declare discrete output metrics, derived from the extrema and temporal means of computed internal states of our model in figure 1.
The contractility,   , and the compliance,   , both control the elastance extrema of the left ventricular and the left atrium.

Measurements
recall, we investigate the identifiability of input parameters assuming certain patient measurements will be available.Here, though, we use of the model solution as a surrogate for the patient data.This allows us to remove issues associated with noise, sampling rate and indeterminacy -we know the answer.Put another way, we consider identifiability of input parameters in the best possible setting, given the available measurements.We utilise the measurements declared in Table 2, which are all clinically plausible.Solving the system, we have acquire time-series waveforms for left ventricular pressure, left ventricular volume, systemic arterial pressure and systemic arterial flow.From these, we compute the derived metrics of left-ventricular stroke volume ( ), e.g.measured by echocardiography [52].Left ventricular pressure can be measured invasively using left heart catheterisation [53], allowing for computation of the pulse pressure in the left ventricle (   ).Systemic arterial pressure and flow are obtained via arterial line measurements and Doppler ultrasound [54,55], leading to evaluation of pulse pressure in the systemic artery (   ) and mean arterial flow ( Q ) respectively.

Computational framework
All model solution and analysis are conducted in Julia [56], which is chosen due to the superior speed and technical versatility in solving differential equation systems [57,58].It has been shown previously that any model analysis involving sensitivity analysis, and by extension identifiability analysis, only becomes too expensive when the model solution time is too large.To solve our single ventricle model we utilise the Vern7 algorithm [59], with each cardiac cycle comprised of 250 time steps and 1 − 6 absolute and relative tolerance.All global sensitivity analysis is conducted with GlobalSensivity.jl[60] and plotting is performed using Makie.jl[61].

Sensitivity analysis
As discussed briefly in Section 2.1, sensitivity analysis is the study of how uncertainty in the output of a model (numerical or otherwise) can be apportioned to different sources of uncertainty in the model input parameters (henceforth termed 'factors') [62].Local sensitivity analysis (LSA) and global sensitivity analysis (GSA) are considered in Sections 3.2.1 and 3.2.3.The Fisher information matrix, utilising the sensitivity matrices, is defined in Section 3.2.2.

Local analysis
Local, derivative based sensitivities are essentially partial derivatives, evaluated at a base state in input parameter space,  0 at time .To compare a parameter 's influence evenly against the output , we scale the raw sensitivity metric by θ ŷ .The result is a relative sensitivity matrix Ŝ with entries Ŝ, .The input parameters  ∈ (1, … , ) (here  = 9) are defined in Table 1 and the measurements  ∈ (1, … , ) (here  = 4) are defined in Table 2.
We define relative sensitivity column vectors associated with a specific model input, , as follows: where Ŝ,1 represents the influence of input parameter  against the measurement 1 (say).To compute the above sensitivity statistics, the inputs are, of course, varied one at a time about  0 .

Fisher information matrix
Another important matrix derived from sensitivity matrices is the square ( × ) Fisher information matrix (FIM) [63]: The FIM is a symmetric matrix representing the information one can extract on factors from the model outputs which correspond to the available measurements [64].It is important to note that the FIM can be constructed from either global or local sensitivity matrices.

Global analysis
We detail each global methods used to quantify parameter influence.Two are variance-based methods (Sobol and eFAST), whereas the Morris method is a one-at-a-time method.To ensure fair comparison between,  = 100, 000 samples are utilised for each method.For the Morris method, this means 100,000 trajectories, taken in the input parameter space.Note that this does not ensure that the same number of model evaluations is needed, but it does ensure equal density of input parameter space is explored for the variance-based methods.As we are interested in the second order indices, given  as the dimensionality of input parameter space (here  = 9), Sobol's method suggest (2 + 1) = 1, 900, 000 evaluations, eFAST requires  = 900, 000 evaluations and Morris' method ( + 1) = 1, 000, 000 evaluations [65,66].To sample the input parameter space bounds of ±15% are prescribed; where a probability distribution is required, a uniform distribution is given with the respective upper and lower bounds.To ensure the true value of the sensitivity indices, we use 'sampling with replacement' -a bootstrapping method (with  = 1000), to calculate confidence intervals for the derived sensitivity indices [67].As Sobol indices are commonly used as a benchmark against other methods, the convergence is only calculated in this case [68,69].
Morris.Morris' method computes the sensitivity of the model outputs, to each parameter   , by sampling the elementary effects (EE hereafter) [70], perturbing one input parameter at a time: , where  is the parameter step size describing the ''levels'' of effects.
Choosing  to be even provides a more symmetric sampling distribution [70], hence we choose  = 100 giving  = 0.51.Morris' method then creates a sampling trajectory, characterised by a single parameter change.One then estimates the mean and variance of the distribution of the |  |, with a high mean (variance) value, implying an important parameter (non-linear effects, or interactions with other inputs).
eFAST.The extended Fourier amplitude test [71] offers a computationally efficient alternative to Sobol indices (below).eFAST utilises mono-dimensional Fourier decomposition along a curve exploring the parameter space.The curve is defined by a set of parametric equations: where   is a transformation function chosen to ensure that the variable is sampled according to the desired probability density function.  is a set of different (angular) frequencies, to be properly selected, associated with each input parameter.As  varies, all the factors change simultaneously along a curve that systematically explores input parameter space.Each   oscillates periodically at the corresponding frequency   , whatever   is.The output  shows different periodicity, combined with the different frequencies   , whatever the model  is.
If the th factor has strong influence on the output, the oscillations of  at frequency   will be of high amplitude.This provides a basis for computing a sensitivity measure, which, for factor   , is founded on the coefficients of the corresponding frequency   and its harmonics.For a full derivation and discussion of frequency choice, see [48,72].
Sobol analysis.Consider discrete outputs,   =  (  ), computed from a sample of those in Eq. ( 1).We regard the   as determined by the choice of input parameters which are continuously distributed (presumably over some physiological range), so one can view the   =  (  ) as integrable multivariate functions.Use a Hoeffding-Sobol decomposition [47], which, for the sake of brevity, we state for three parameters: If we assume that input parameters are independent and can be transformed to   ∈ [0, 1] are normalised and that the eigenfunctions in the expansion on the right hand side exhibit the orthogonality property: then it is straightforward to show, by a recursive process of integratingout first all, then two, then a single variable, that these eigenfunctions have a simple and intuitive realisation: Above, (..) denotes the expectation value and , ,  = 1, … , 3 and    represents the complement set of the input parameter   .Taking a relative variance yields a single intuitive measure-the first order Sobol index of the parameter   : where  is the expectation operator.The inner expectation operator functions such that the mean of  is taken over all possible values of    while keeping   fixed.The outer variance is taken over all possible values of   .Then utilising the identity [73]: where Var   (    ( |  )) measures the first order (additive) effects of   on the model outputs.Another popular variance measure (prominent in this work) are the total order estimators, first introduced by Homma [74]: Here   , measures the total effect, i.e., first and higher order effects )) must give the contribution of all terms in the variance decomposition, which include the input   .
To interpret Sobol indices, we note that a quantity of interest with high values of either index, for a given factor, suggests that measurement of that quantity may provide substantial information about that parameter.If   is low but   , is large, then parameter   impacts the quantity of interest primarily through interactions with other parameters and so still may substantially affect the quantity of interest.

Average parameter influence
So far, we have not postulated a metric for the overall influence of an input parameter across all the outputs.In [45], Li et al. derive such a metric based upon the FIM  in Eq. ( 9).Li's method has only been applied to an FIM derived from LSA.We now extend to an FIM derived from GSA, deriving a parameter influence, or effect using principal component analysis (PCA) [45,75].The principal components (PC) are the eigenvectors of the FIM.
Let  be the matrix of the ordered PC eigenvectors of  , in which the absolute value of each element   reflects the contribution of the th parameter to the variance of the th output.We follow Li et al. [45], who measure an overall effect for the th parameter as: where 0 ≤   ≤ 1 and   represents the non-zero eigenvalues of  .This measure reflects the difficulty in determining the th factor when only a single factor is estimated.Accurate estimation from limited data sets is favoured by large values of   .

Orthogonality analysis
Using Eq. ( 8), we can define the orthogonality between two input parameters by utilising the derived sensitivity vectors.  can be viewed as the averaged orthogonality between two input parameters   and   , aggregated across all the outputs (see Table 2).
where, ‖.‖ denotes an Euclidean norm, the sin function ensures orthogonality measure   ∈ [0, 1].This measure of orthogonality will be utilised below, for our extended subset selection method.It will also be used to rank parameters based on the most independent effects while not considering influence at all.To calculate a rank based on orthogonality, we take the mean orthogonality score for an input parameter with respect to all the others across all the outputs.

Parameter identifiability
We require a practical strategy for finding suitable input parameter subsets for the purpose of model personalisation.Our method is based upon a technique dating from 2004, due to Li et al. [45], originally applied to local sensitivity of bioreactor design; it offers an intuitive balance of parameter influence and orthogonality.This selection will use an aggregated identifiability index, which is a simple, scalar product of a measure for: (i) parameter influence given in Eq. ( 17) and (ii) parameter linear independence given in Eq. ( 18) (we want a subset, the members of which are optimally linearly independent): Above,   is our measure of influence obtained from Eq. ( 17) and   is our measure of the linear independence obtained from Eq. ( 18).
The concept of orthogonality underlies the method of Li et al. [45], as follows.Parameter dependence is quantified from the FIM (Eq.( 9)).The rank of  , defined as the dimension of the vector space spanned by its columns [76], gives the number of identifiable combinations of input parameters at any given model operating point [77,78].Given that we have 4 measurements, 4 identifiable (independent) parameters are required to span the whole output space.It is thus our extension to perform global sensitivities to identify which input parameters may be uniquely extracted from the available measurements.
Assuming ŝ , for  = 1, … ,  with  <  are all linearly independent, we find the projection of another vector, ŝ, into the subspace spanned by ŝ where  < , effectively by removing out its orthogonal projection (which lies outside that subspace).The remaining part is given by Consider, now, a new candidate sensitivity vector ŝ , for a generic th input parameter.The extent to which ŝ is linearly dependent upon the already-chosen ŝ , is measured by finding the above projection of ŝ , onto the subspace spanned by the ŝ , that is, by removing from ŝ its orthogonal compliment, ŝ⟂  .Accordingly, ŝ∥  is defined by its expansion coefficients,   .The latter may, in fact, be efficiently computed in an optimisation process, with solution [45]: We summarise our algorithm for input parameter subset selection.
1.For each parameter,   ,  = 1, … , , each having relative sensitivity vector ŝ , calculate its overall effect, using Eq. ( 17); 2. Select the parameter with the highest value of   ,  = 1, … , , to be the first parameter in the selected set; 3.For  < , repeat the following steps until no more parameters can be added to the accumulating set.For the th candidate: (a) Use Eq. ( 21) to find the nearest vector ŝ∥  , to the present candidate, lying in the subspace already spanned by the  (say) currently selected parameters.(b) Use Eq. (18) to calculate the orthogonality between ŝ∥  and ŝ as follows (Of course,   is a proxy for the magnitude of the orthogonal projection of ŝ which, in turn, measures an overall orthogonality for candidate  relative to the  already-selected parameters.)(c) Attribute to candidate parameter  a simple aggregate identifiability index which reflects both its sensitivity and orthogonality (d) Provided   > 0.05 [32], include in the set that parameter with max  (  ).
4. If  ≥  form all (m-1) tuples of the already selected parameters.The number of possible candidates is Use Eq. ( 22) to calculate the orthogonality of the input parameter   across all  possible combinations of parameters  , .Determine the worst case scenario (  = min( , )) and continue with the calculation of   .
5. Continue until no more parameters (elements) can be added.
H. Saxton et al.

Hyperspace dimension
From Eq. ( 16), one can interpret the total order Sobol indices as [79]: i.e., for a given input parameter   , the total order indices are the first order indices plus all higher order interactions.Subtracting the first order indices where we denote  , as the higher order interactions, for an input parameter   .If  , ≈ 0.0, one can deduce that the model inputs impact the model outputs in an independent way.Total order sensitivity indices have been shown to exhibit superior ability to recover the true value of sensitivity indices, when compared to first order indices [69,80].Thus, one can utilise the total order indices to quantify the complexity of the input parameter space.We now define the investigative test procedure for our problem as follows: 1. Perform an extensive SA at the maximum bounds for our study: ±15% to compute  1 ,  2 and   .Calculate the respective confidence intervals.2. Ensure all sensitivity indices exhibit a statistical variation of less than 5% of the mean value.3. Calculate   .
(a) if   < 0.01, proceed using only   .(b) If   > 0.01, proceed as follows but examine  1 and   to ensure each sensitivity index has converged.
4. Ensure consistent sampling density and converged sensitivity indices, iteratively reduce the hyperspace dimension, recording the rank and sensitivity values of each input parameter.5. Once the hypercube dimension investigated is less than ±0.05% from base state, finish the investigation.6.Once all input parameter rankings and values have been recorded, examine the variation between the edges of the input parameter hypercube compared to the minimal variation from the base state.
If a consistent ranking is exhibited at all hypercube sizes, one can infer a less complex input parameter space, with obvious consequences model personalisability.On the other hand, if we observe large variations in input parameter rankings when the hypercube sizes are varied, one would infer a complex input parameter space hyper-surface -a greater encumbrance to personalisation (due to multiple possible local minima causing the input parameter ranking variations), when compared to the previous example.

Results
Sections 4.1 and 4.2 respectively examine LSA and GSA results for our model and Sections 4.2.1-4.2.3 the Morris, eFAST and Sobol GSA methods explored.In Section 4.3 the overall effect of input parameters is compared between methods utilising Eq. ( 17).We then compare the ranking of input parameters, based solely on their orthogonality score, between methods utilising Eq. ( 22), in Section 4.4.Section 4.5 examines the stability of input parameter identifiability when we extend subset selection to GSA methods and finally Section 4.6 examines the complexity of the input parameter space of our single ventricle model, using the method of Section 3.6.

Local sensitivity
Fig. 2 is the local sensitivity matrix.The minimal elastance of the left ventricle   is most influential, across all measurements, followed by the maximal elastance   and windkessel factors   ,   and   .The end diastole timing parameter   and the valve parameters   and   are the least influential.

Global sensitivity
The mean vs. variance plots are displayed, for Morris' method in Fig. 3. Figs. 4 and 5 show first order indices (panel A), total order indices (panel B) and the higher order indices (panel C) as defined in Eq. ( 24), for the eFAST and Sobol methods, respectively.

Morris method
Fig. 3, plotted on a log 10 scale, shows using Morris' method, that both valve parameters   ,   and the minimal elastance   have high mean and variance values against all 4 measurements, implying these inputs are influential and have either a non-linear relationship with the output, or non-linear interactions with other inputs.Here the venous compliance   has a low mean and variance for all 4 measurements, implying   has little influence on the output, and may be fixed.

eFast method
The sensitivity indices generated from the eFAST method are displayed in Fig. 4. Panel C displays the higher order sensitivity indices  ℎ ≈ 0 with the highest order interaction value of 0.0066, for   , impacting the pulse pressure of the systemic artery.Thus, using the eFAST method, we infer the inputs to act independently on the outputs.Examining panels A and B, the ventricular elastance parameters   ,   appear influential across all 4 measurements.The systemic resistance   appears influential across all measurements apart from the ventricular pulse pressure.The arterial compliance   appears influential when the measurement includes a pressure.The other system parameters appear to have little influence across all outputs.

Sobol indices
Fig. 5 displays the Sobol indices.As in the eFAST case, we note that panel C, for the higher order indices, are all approximately 0. The largest higher order interaction is the minimal elastance   , with a value  ℎ = 0.0099, impacting the stroke volume of the left ventricle.Because the higher order indices are of very low value, the second order indices are not displayed here but in Fig. 7 in the appendix instead.Panels A and B show the first and total order indices (their respective convergences are displayed in Figs. 8 and 9).  appears most influential across all measurements, with the system parameters   ,   and   next.

Input parameter influence comparisons
Figs. 2, 3, 4 and 5 show our sensitivity matrices for each input parameter, using different methods.Utilising the method in Section 4.3, Table 3 displays the average influence ranking of all input parameters, averaged across all 4 measurements.Although exact influence values differ, all sensitivity measures rank the minimal elastance   as the most influential across all measurements.All sensitivity measures, except the Morris' method, rank the arterial compliance   as second most influential, with Morris' method attributing an influence measure and order of magnitude lower than all other methods.Interestingly, the Morris sensitivity measure ranks the valve parameters as the next most influential.All global measures apart from Morris' rank input parameters in the same orders, apart from the parameters with negligible influence values.The first order and total order indices exhibit the same ranking, which once again is indicative of a system driven by independent input factors  ℎ ≈ 0. The Local sensitivity matrix displays a very similar ranking to the global measures, with small differences in positions.

Orthogonality analysis
Fig. 6 (Panels A-F) display the orthogonality matrices for the respective sensitivity measures.Overall, the orthogonality rankings of input parameters are different using different sensitivity measures.In some instances, for example, for the eFAST and Sobol methods, there are vastly different orthogonality scores between the first and total order indices, e.g.  .Despite this, there are still common themes for the orthogonality scores:   being consistently independent,   having independent effects from   , and   ,   having dependent effects.Panels G-L underscore varied orthogonality scores between different sensitivity measures.Sobol first order, Sobol total order, eFAST first order and the local measure are consistent with more orthogonal parameters.Morris' and eFAST total order produce orthogonality rankings which suggest that the input parameters are more dependent on each other.
Table 4 displays the average (as defined in Section 3.4) input parameter ranking based on orthogonality.In contradistinction to the influence case shown in Table 3, no clear patterns in the ranking emerge.  appears the most orthogonal for the eFAST first order, Sobol first order and Sobol total order methods.The venous compliance   ranks as least orthogonal in all sensitivity measures, apart from the eFAST total order and the Sobol first order, but even in these settings, the rank of   is low.Examining the parameter orthogonality value column in Table 4, we see a large variation in the average values of orthogonality for each input parameter, with the lowest ranked parameter exhibiting an average orthogonality score between (0.21-0.549).

Identifiability analysis
Utilising Eq. ( 9), for every sensitivity matrix generated from LSA and GSA, we find the FIM to be singular, i.e., certain model parameters' effects are totally dependent on others.( ) = 4 for the single ventricle model defined in Eqs. ( 2)-( 6) implies 9−4 = 5 non-identifiable parameter combinations.Thus, we expect 4 identifiable input parameters.From the definition of the identifiability index (Eq.( 19)), we use a cut off of  < 0.05 [32] although further investigation of the reliability of this measure is indicated.We note all sensitivity methods, except Morris' obtain the expected number of identifiable input parameters in Table 5.All global methods identify the minimal elastance (  ) and arterial compliance (  ) as the most identifiable parameters.The local sensitivity measure agrees with the global methods, with the minimal elastance as the most identifiable and the arterial compliance as identifiable (but this is ranked 4th).Sobol first order indices and the local measure are the only methods that find the systemic resistance   as identifiable.eFAST and the Morris method first order indices are the only methods to find the maximal elastance as also identifiable.Both total order indices for the global methods give the same set of identifiable input parameters.Examining the identifiability index value of each input parameter, all global methods except Morris' attribute a similar value for each rank position, indicating that although the input parameter rankings may vary between methods, their quantifiable identifiability value remains constant.Fig. 6.Orthogonality Matrix and Orthogonality histograms: Panels A-F show the orthogonality matrices for the local, Morris, eFAST first order, eFAST total order, Sobol first order and Sobol total order methods, respectively.A value of 1(0) indicates that the two input parameters have orthogonal effects on across all the outputs (contribute the same effect on the output).Panels G-L are histograms of the respective orthogonality matrices, indicating the distribution of orthogonality present within the input parameters when computed through the different sensitivity measures.

Table 4
Parameter orthogonality ranking: A table displaying the rank of input parameters based on their average orthogonality score, calculated by taking the mean orthogonality score for each input parameter across all outputs for each sensitivity measure.

Sensitivity metric
Parameter

Hypercube dimension
Next, in Table 6, extending the input parameter space volume, from local to one characterised by a parameter variation of ±15% reveals an input parameter ranking which remains constant from a boundary ±0.01%upward.When using the local sensitivity measure, we note that the ranking is similar to the global setting, with the minimal elastance   and arterial compliance   ranking first and second.The other input parameters, when using the local measure, appear to vary by only a single rank position, when compared to the global setting.As each hypercube dimension is sampled with the same density, one would expect to see the parameter influence value to remain the same.Here, as the hypercube dimension is extended, we observe some slight variation in the influence value.This shows that, as the hypercube is extended, there is some quantifiable change in an input parameter's influence.

Discussion
With an aim to address the stability of input parameter identifiability, we extended the subset selection method of Li et al. [45] to global sensitivities, then utilised different global sensitivity methods to interrogate the input parameter space.Tables 7 and 8 highlight issues surrounding this problem.Comparing influence and orthogonality rankings directly, we noted that influence had a much more consistent ranking with all methods, except the local and the Morris methods.When the ranks are based on orthogonality, we observe that no methods exhibit a consistent ranking of input parameters.This is further apparent when the extended parameter subset selection method is applied.Table 5 highlights that all global methods except Morris' found the minimal elastance   , the arterial compliance   and the venous compliance   to be identifiable, however the 4th parameter found to be identifiable varied between each method.The total order indices for eFAST and Sobol produced the same subset of identifiable parameters in the same order, which is reassuring, although we see that the model is driven mainly through independent effects ( ℎ ≈ 0), where the total order indices capture all contributing affects to the output variance.Thus, the identifiable subset returned is the same.
It is clear from the way each sensitivity method aggregates that orthogonality has a large impact on an input parameter's identifiability and therefore should be examined further.However, it is important to note that the extended subset methodology (with results in Table 5) utilises the concept of orthogonality differently to the way it was analysed in Table 4. Once the number of the selected input parameters is greater than the number of measurements available, the orthogonality score used to calculate the identifiability index from Eq. ( 19) (i.e., the worst case/maximum orthogonality) is chosen.Thus, for input parameters deemed unidentifiable, the rank should not be examined too closely, as there are an infinite number of unidentifiable ranking positions for the input parameters.Moreover, with the extended subset methodology, we examine orthogonality against groups of input parameters, so the rankings presented in Table 4 (which are based on averages of orthogonality pairings of the parameters) may not translate directly to the extended subset methodology.
The Morris method has failed to return consistent results for both the influence and an identifiable subset.While it is still popular for higher dimensional models, its low ability to explore input parameter space has here been highlighted.This ability deteriorates exponentially with increasing dimensions [81].Ideally, one would use a variance based method such as eFAST or Sobol indices to characterise input parameter, however this is often not utilised due to the associated computational expense.Sensitivity analysis is driven by the speed in which a dynamical system can be solved.For our system, see Section 3.2.3.The time taken to compute Morris, eFAST and Sobol indices on 28 threads was 3568.6, 3325.1 and 6893.0 s respectively.Thus, eFAST presents itself as a reliable and efficient GSA method.Note, for the Sobol method, the time quoted above included the computation of second order indices also; if we were only interested in the first and total order the number of model evaluations would be the same as Morris' method.Thus, given one can optimise the model solution time, efficient GSA is assured.This may be achieved through use of surrogates [82,83] or by utilising the efficient ODE solvers in DifferentialEquations.jl[58], as here.Overall, eFAST is a reliable method to assess the uncertainty of the single ventricle model.However it relies on a sinusoidal function to sample input parameter space, which creates two problems.This sampling method produces a zigzag pattern in the input parameter space and it can struggle to capture the extremes of the input parameter space [84].When the size of the input parameter space increases, and with more physiologically detailed models, the eFAST method may also struggle to return true input parameter sensitivities.On the other hand, the Sobol methodology utilises Quasi-Monte Carlo sampling strategies, which allow for an easy computation of the confidence intervals associated with the sensitivity index.No such method exists for the eFAST methodology, due to the sampling nature of the method [85].Therefore, for assessing uncertainty, the Sobol method is still preferred, because of its ability to compute confidence intervals, through the bootstrapping methodology alongside the sensitivity indices.
A secondary aim of our investigation was to develop a GSA based methodology for a mapping of the input parameter space.We have shown in Table 6 that as one migrates from base state, input parameters' influence rank remain constant, when exploring an extended region.A possible explanation for the change in parameter rank moving from base state to ±0.01% is that local SA is a linearisation around a point meaning that any non-linear effects in the model are neglected.The total order Sobol indices method captures all effects -no matter how minor.In addition, in Table 6 the value quantifying parameters' influence remains largely constant with minor variations, implying that as the size of the hyperspace is extended, new domains of parameter influence are reached.However, the effects as a whole do not change because of the extension of the hyperspace dimensions.This indicates a largely additive model, representing a flat input parameter hyperspace.A (mostly) additive nature suggests that the single ventricle model is a good candidate for personalisation.Validation of this conclusion can be performed from patho-physiological patient measurements; however, our present model is only shown to be effective in identifying patient's characteristics from exercise data using local optimisation [86].Our findings, utilising the hyperspace dimension test are consistent with this.
Despite the promising findings of our hyperspace dimension test, there are important issues to consider.First, we use a GSA method is used to quantify the input parameter space.Often, GSA is performed without any assessment of the error associated with the sensitivity value.Previous work has shown that in order to achieve the true input parameter influence, convergence must be achieved.Thus, we propose Sobol indices are used for the hyperspace dimension test.Second, one should ensure the same density of sampling is applied at every hyperspace dimension size.The selection of samples suffers 'the curse of dimensionality' [87].But if a large enough sample size over the hypercube was not utilised, consistent sample sizes close to the base state would become too small.Here, we have utilised the total order indices due to  ℎ ≈ 0. If this was not the case, one would have to perform investigations with both the first and total order indices.The first order indices are much harder to converge [69,79,80], thus for a model with slightly more complexity and higher non-linearity, it may not be possible to perform such a convergence test.However, for dynamical systems where the outputs are driven by independent effects of the inputs, this test will prove useful.
For the personalisation process, investigations must be assumed to rely on a fixed set of measurements.In reality, one may not have access to a rich and diverse set.While our methodology defines the ideal offline scenario to investigate identifiability -before using live patient data -it should be acknowledged ever investigation is constrained by the measurements available.So, every time a new measurement becomes available, a new identifiability study is indicated.This also has implications for parameter inference, because parameter identifiability is only valid in the presence of a set of measurements.

Conclusion
Our study provides a clear and intuitive investigation of a key stage of the personalisation process.We have studied a single ventricle, 9 parameter, 0D model, to probe the identifiability of its input parameters, in the presence of 4 synthetic clinically available measurements.We have: (i) extended the parameter subset selection method of Li et al. [45], to encompass the global nature of the personalisation process, (ii) revealed how a new and different set of globally identifiable input parameters could be obtained and (iii) provided novel perspectives, compared to the previous local studies.Assessing the stability of this identifiable input parameter subset, we employed various global and local measures of input parameter sensitivity, revealing how alternative sensitivity methods which depict input parameter space in contrasting ways lead to similar but subtly different identifiable input parameter subsets (driven mainly by the dissimilar orthogonality between input parameters).Finally, we have detailed a novel and intuitive input parameter hypersurface structure investigation, utilising Sobol indices.The connection with Sobol index error evaluation provides a guide for mapping of the complexity of input parameter space, with a view to aid the inverse problem.When applied to the single ventricle model, within the presence of the 4 chosen measurements, the single

Fig. 2 .
Fig. 2. Local relative sensitivity Matrix: Shows the local relative sensitivity matrix, measuring input parameters' (column headings) influence on specific model outputs (row headings).

Fig. 3 .
Fig. 3. Morris' method scatter plots: Each plot displays a normalised mean value plotted against the variance value for each input parameter on a log 10 scale.Panel A: Morris' method results for the stroke volume of the left ventricle.Panel B: Morris' method results for the pulse pressure.Panel C: Morris' method result for the pulse pressure in the left ventricle.Panel D: Morris' method results for the mean systemic flow.

Fig. 4 .
Fig. 4. eFAST sensitivity matrices: Each matrix, with input parameters as column headings and specific model outputs as row headings, displays an influence value for an input parameter against a specific output.Panel A: the first order indices.Panel B: the total order indices.Panel C: the difference sensitivity matrix as defined in Eq. (24).

Fig. 5 .
Fig. 5. Sobol sensitivity matrices: Each matrix, with input parameters as column headings and specific model outputs as row headings, displays an influence value for an input parameter against a specific output.Panel A: the first order indices.Panel B: the total order indices.Panel C: the difference sensitivity matrix, as defined in Eq. (24).

Fig. 7 .
Fig. 7. Second Order Sobol indices: The second order Sobol indices are presented as lower triangular matrices due to their symmetric nature.Each matrix element displays the influence an interaction between two input parameters have on a selected output.Panel A: Displays the second order indices stroke volume for the left ventricle.Panel B: Displays second order indices for the pulse pressure for the left ventricle.Panel C: Displays the second order indices pulse pressure for the systemic artery.Panel D: Displays the second order indices for the mean systemic flow.

Quantification of Input Parameter Space Complexity: We
detail an investigative test, driven by global sensitivity analysis, to quantify the complexity of the input parameter space.

Table 1 Input parameters for the single ventricle model.
Each input parameter's unit is stated alongside a chosen initial value for the 9 parameter, single ventricle model. is the cardiac cycle length and is fixed such that  = 1 s.The ventricular shift parameter  shift = 0 s as no atrium is present.
multiplicative interactions) of input parameter   .One can consider this by recognising that Var    (   ( |   )) is the first order effect of    , so Var( ) − Var    (   ( |

Table 3
Parameter influence ranking: A table displaying the ranking of each input parameter influence, averaged across all 4 measurements.Rankings are displayed for both local, global, first order and total order sensitivity measures.

Table 5
Identifiable input parameters: Table displaying the identifiable input parameters calculated using the global subset selection method.Parameters in red indicate an unidentifiable input parameter utilising a cut off of  < 0.05.

Table 7
[45]mean rank and the range of the input parameters: Table showing the effect of different parameter subset methodologies (influence, orthogonality and the extended Li methodology[45]) when we stratify across all sensitivity metrics.

Table 8
The mean rank and the range of input parameters across all subset selection methodologies when we stratify by different sensitivity methods.