Sensitivity Analysis: A Useful Tool for Bifurcation Analysis

Abstract

Sensitivity analysis and bifurcation analysis are closely related to each other. In sensitivity analysis, especially global sensitivity analysis the effects of input parameter spaces on output quantities of interest are studied. On the other hand, in bifurcation analysis the critical points within feasible regions of parameters are detected where the long-term dynamics changes qualitatively. Prior to bifurcation analysis, it is important to identify the bifurcation parameters. In complex and computationally expensive problems which consist plenty of uncertain parameters, it is essential to find a set of bifurcation parameters before bifurcation analysis. Global sensitivity analysis is a powerful tool to identify the bifurcation parameters which contribute most on output uncertainty. Global sensitivity analysis is the first step toward bifurcation analysis which helps in dimension reduction during bifurcation analysis. As an example, in this chapter, a multi compartment, lumped-parameter model of an arm artery is considered and global sensitivity analysis (Sobol’s method) is applied to identify the bifurcation parameters of the arm arteries.

1. Introduction

Sensitivity analysis and bifurcation analysis are closely related to each other. In sensitivity analysis, we study how the uncertainty in the output of a mathematical model or system (numerical or otherwise) can be apportioned to different sources of uncertainty in its inputs [1]. On the other hand, in bifurcation analysis the critical points within the feasible regions of parameters are detected where the long-term dynamics changes qualitatively [2]. Prior to the bifurcation analysis, it is important to identify the bifurcation parameters in complex and computationally expensive problems that consist plenty of uncertain parameters. Sensitivity analysis is a powerful tool to identify the bifurcation parameters which contribute most on output uncertainty. Also, sensitivity analysis helps in dimension reduction during the bifurcation analysis by fixing less influential parameters on their nominal values.

Sensitivity analysis can be divided into two categories, local sensitivity analysis (LSA) and global sensitivity analysis (GSA). In LSA a parameter value is perturbed around its nominal values at a time, keeping other parameters fixed on their nominal values. The procedure is repeated for all parameters one by one to study their impact on output variables. LSA techniques are simple, easy to implement and computationally less expensive. On the other hand, LSA is not suitable for non-linear models and does not explore the impact of entire parameter spaces and their interactions effects on output variables [3, 4]. In order to overcome the limitations associated with the LSA, GSA can be used. In GSA, the analysis is performed over entire feasible regions of the input parameters and quantifies the impact of parameter interactions on output variables. The only deficiency related to the GSA is its computational cost [5, 6, 7, 8, 9, 10, 11, 12] (Figure 1).

Figure 1.

A simplified 5-step procedure to identify the bifurcation parameters using global sensitivity analysis.

Figure 2.

Model geometry of arm artery (A), with total number of arterial segments, Ns = 15, including 12 non-terminal and 3 terminal segments. Each non-terminal and terminal segment is represented by its corresponding non-terminal (C) and terminal electrical circuits (D). Pressure waves is used as an input boundary condition (B) and pout = 15 mmHg which is mean venous pressure used to calculate boundary outflow. The parameter values of each arterial segment are given in Table 1.

In this chapter, the main questions of interest are:

1. How to identify the bifurcation parameters in a model having plenty of input parameters?

2. Which parameters could be exempted from the bifurcation analysis (dimension reduction)?

This chapter seeks to answers the above-mentioned questions using a simplified 5-steps procedure of uncertainty and sensitivity analysis. As an example, a multi-compartment, lumped-parameter model of arm arteries is considered [4] and global sensitivity analysis (Sobol’s method) is applied to identify the bifurcation parameters (electrical) of the arm arteries.

2. Lumped-parameter model of the arm arteries

In this section, the major arteries of the arm are divided in to number of non-terminal and terminal arterial segments (nodes). The total number of arterial segments, ${\mathrm{N}}_s=15$ including 12 non-terminal and 3 terminal segments. Each non-terminal and terminal arterial segment is represented by its corresponding non-terminal and terminal electrical circuit.

Applying Kirchhoff’s current and voltage laws on electrical representation of arm arteries, the following mathematical equations for pressure and flow are obtained:

Pressure and flow equations at non-terminal nodes:

Flow equation:

$$\begin{array}{l}{\dot{q}}_i=\frac{p_{i-1}-p_i-R_i{q}_i}{L_i},i=1,2,3,\dots 15\text{and}i\ne 11,13\\ {\dot{q}}_{11}=\frac{p_6-p_{11}-R_{11}{q}_{11}}{L_{11}}\\ {\dot{q}}_{13}=\frac{p_{11}-p_{13}-R_{13}{q}_{13}}{L_{13}}\end{array}$$ E1

Pressure equation:

$$\begin{array}{l}{\dot{p}}_{\mathit{i}}=\frac{q_i-q_{i+1}}{C_i},i=1,2,3,\dots 15\text{and}i\ne 6,11\\ {\dot{p}}_6=\frac{q_6-q_{11}-q_7}{C_6},\hspace{1em}{\dot{p}}_{11}=\frac{q_{11}-q_{12}-q_{13}}{C_{11}}\left(\mathrm{at}\text{bifurcation}\right)\end{array}$$ E2

Pressure and flow equations at terminal nodes:

$$\begin{array}{l}{\dot{q}}_{\mathit{in}}=\frac{2p_{\mathit{\text{in}}}-2p_i-R_i{q}_{\mathit{\text{in}}}}{L_i}\\ {\dot{p}}_i=\frac{q_{\mathit{\text{in}}}-q_{\mathit{\text{out}}}}{C_i}\\ {\dot{q}}_{\mathit{\text{out}}}=\frac{2p_i-2p_{\mathit{\text{out}}}-2R_b{q}_{\mathit{\text{out}}}}{L_i},i=10,12,15\end{array}$$ E3

where, $$R_i,C_i\text{and}{L}_i$$ is the blood flow resistance, compliance of the vessel and blood inertia of $$i^{\mathit{th}}$$ segment of the arm arteries respectively. The electrical parameters $$\left(R_i,C_i,L_i\right)$$ of $$i^{\mathit{th}}$$ segments are related with structural parameters $$\left(E_i,l_i,d_i,h_i\right)$$ as,

$$R_i=\frac{8\nu l_i}{\pi {\left(\frac{d}{2}\right)}^4},C_i=\frac{\rho l_i}{\pi {\left(\frac{d}{2}\right)}^2},L_i=\frac{2\pi {\left(\frac{d}{2}\right)}^2{l}_i}{E_i{l}_i}$$ E4

where, $$E_i$$ is the Young modulus, $$l_i$$ denotes length of the vessel, $$d$$ is the diameter of the vessel and $$h_i$$ represents the wall thickness of $$i^{\mathit{th}}$$ segment of the vessel. Moreover, $$\nu$$ (0.004 Pa s) is the blood viscosity and $$\rho$$ ( $$1050 {\mathit{\text{kgm}}}^{-3}$$ ) is the blood density. The nominal values of all parameters of arm segments are given in Table 1. The geometry along with the values of the parameters is taken from [13, 14].

Nodes	E	l	d	h	R	C	L
units	$${\mathit{kgm}}^{-2}{s}^{-2}$$ $$\times {10}^5$$	$$m$$ $$\times {10}^{-2}$$	$$m$$ $$\times {10}^{-3}$$	$$m$$ $$\times {10}^{-4}$$	$${\mathit{kgs}}^{-1}{m}^{-4}\times {10}^6$$	$${\mathit{kg}}^{-1}{s}^2{m}^4\times {10}^{-11}$$	$${\mathit{kgm}}^{-4}\times {10}^6$$
1	4	6.1	7.28	6.2	3.539	7.454	1.539
2	4	5.6	6.28	5.7	5.868	4.778	1.898
3	4	6.3	5.64	5.5	10.15	4.035	2.648
4	4	6.3	5.32	5.3	12.82	3.514	2.976
5	4	6.3	5	5.2	16.43	2.974	3.369
6	4	4.6	4.72	5	15.10	1.9	2.76
7	8	7.1	3.48	4.4	78.90	0.667	7.838
8	8	7.1	3.24	4.3	105	0.531	9.042
9	8	7.1	3	4.2	142.9	0.448	10.55
10	8	2	2.84	4.1	55.11	0.1207	3.647
11	8	2	4.3	4.9	31.94	1.067	4.844
12	16	6.7	1.82	2.8	1173	0.0834	31.88
13	8	7.9	4.06	4.7	40.19	0.9366	5.434
14	8	6.7	3.48	4.6	50.22	0.80	6.075
15	8	3.7	3.66	4.5	33.60	0.3958	3.693

Table 1.

Numerical values of parameters for each node of the arm arteries (shown in Figure 2).

The value of boundary resistance $$\left(R_b\right)$$ on three terminal nodes is $$3.24\times {10}^9 {\mathit{kgs}}^{-1}{m}^{-4}$$ , $$\nu =0.004 {\mathit{kgs}}^{-1}{m}^{-1}$$ and $$\rho =1050 {\mathit{kgm}}^{-3}$$ [4, 5, 13, 14].

3. Uncertainty and sensitivity analysis

Uncertainty analysis (UA) and sensitivity analysis (SA) are closely related; however they represent two different disciplines. Uncertainty analysis assesses the uncertainty in model outputs caused by uncertainty of its inputs. Whereas, sensitivity analysis study the impact of input quantities of interest (QoI) on output quantities of interest (QoI). In this study, the input (QoI) are electrical parameters ( $$R_i,C_i,L_i$$ ) and output (QoI) are pressure and flow at each node of the arm arteries. Further, for uncertainty analysis Latin hypercube sampling (LHS) is used and variance-decomposition method (Sobol’s method) is used for global sensitivity analysis (GSA).

Compared to the high-dimensional cardiovascular models (3D, 2D, 1D), lumped-parameter models of the cardiovascular system (CVS) are computationally less expensive, therefore they are suitable for GSA. In our previous studies, we found that for lumped-parameter models of the CVS, the Sobol’s method is computationally less expensive as compared to the other variance-decomposition methods, like FAST and sparse grid stochastic collocation method based on Smolyak algorithm [5, 12].

3.1. The method of Sobol

The method of Sobol is the variance-decomposition method used for global sensitivity analysis. The method decomposes the output variance of a system or model into fractions and assigns them to the inputs factors. For example, given a model of the form $$Y=f\left(X\right)=f\left(x_1,x_2,\dots ,x_k\right)$$ , where $$X$$ is the vector of $$K$$ uncertain parameters, which are independently generated within a unit hypercube i.e. $$x_i\in {\left[0,1\right]}^k$$ for $$i=1,2,3,\dots ,K$$ . Compared to the other GSA methods, the Sobol’s method is one of the most commonly used variance-decomposition method, because of its ease of implementation. The method is primarily based on the decomposition of output $$Y$$ into summands of elementary functions in terms of increasing dimensionality [1, 8],

$$f\left(x_1,x_2,\dots ,x_k\right)=f_0+\sum _i^k{f}_i\left(x_i\right)+\sum _i^k\sum _{i<j}^k{f}_{\mathit{ij}}\left(x_i,x_j\right)+\dots +f_{1,2,3,\dots ,k}\left(x_1,x_2,x_3,\dots ,x_k\right)$$ E5

In Eq. (5), $$f$$ is integrable, $$f_0$$ is a constant, $$f_i$$ is a function of $$x_i$$ , $$f_{\mathit{ij}}$$ is a function of $$x_i$$ and $$x_j$$ and so on. Furthermore, all the terms in the functional decomposition are orthogonal, which leads toward the following definitions of the terms of the functional decomposition in term of conditional expected values.

$$\begin{array}{l}{f}_0=E\left(Y\right)\\ f_i\left(x_i\right)=E_{x\sim i}\left(Y|x_i\right)-f_0\\ f_{\mathit{ij}}\left(x_i,x_j\right)=E_{x\sim \mathit{ij}}\left(Y|x_i,x_j\right)-f_0-f_i-f_j\\ \dots \end{array}$$ E6

where, $$E$$ describes the mathematical expectation and $$x_{\sim i}$$ denotes all parameters except $$x_i$$ and so on. The total unconditional variance can be obtained by,

$$V=\underset{{\mathrm{\Omega }}^K}{\int }{f}^2\left(X\right)\mathit{dx}-f_0^2$$ E7

From Eq. (7), the total unconditional variance can be decomposed in a similar manner like in Eq. (5) as,

$$V=\sum _i^k{V}_i\left(x_i\right)+\sum _i^k\sum _{i<j}^k{V}_{\mathit{ij}}\left(x_i,x_j\right)+\dots +f_{1,2,3,\dots ,k}\left(x_1,x_2,x_3,\dots ,x_k\right)$$ E8

where, $$V$$ is the variance operator. The relationship between functions and partial variance are given by,

$$\begin{array}{l}{V}_i=V_{x_i}\left(E_{x\sim i}\left(Y|x_i\right)\right)=V\left(f_i\left(x_i\right)\right)\\ V_{\mathit{ij}}=V_{x_i,x_j}\left(E_{x\sim \mathit{ij}}\left(Y|x_i,x_j\right)\right)-V_i-V=V\left(f_{\mathit{ij}}\left(x_i,x_j\right)\right)\\ \dots \end{array}$$ E9

Dividing both sides of the Eq. (8) by $$V$$ , we get:

$$1=\sum _i^k{S}_i\left(x_i\right)+\sum _i^k\sum _{i<j}^k{S}_{\mathit{ij}}\left(x_i,x_j\right)+\dots +S_{1,2,3,\dots ,K}\left(x_1,x_2,x_3,\dots ,x_K\right)$$ E10

Where,

$$\begin{array}{l}{S}_i=\frac{V_{\mathit{ij}}}{V}, \text{and}\\ S_{\mathit{ij}}=\frac{V_{\mathit{ij}}}{V}\end{array}$$ E11

where, $$S_i$$ is the main effect (first order sensitivity index) of the $$i^{\mathit{th}}$$ parameter on output uncertainty and $$S_{\mathit{ij}}$$ is the interaction effect of $$i^{\mathit{th}}$$ and $$j^{\mathit{th}}$$ parameters on output uncertainty. Further, the total sensitivity index, $$S_{T_i}$$ can be calculated as,

$$S_{T_i}=\frac{E_{x\sim i}\left(V_{x_i}\left(Y|x\sim i\right)\right)}{V}=1-\frac{V_{x\sim i}\left(E_{x_i}\left(Y|x\sim i\right)\right)}{V}$$ E12

In general, the main effect is used identify the most influential parameters (bifurcation parameters) and the total effect is taken into account for those parameters which are exempted from bifurcation analysis (factor fixing). The total effect, $$S_{T_i}$$ of the $$i^{\mathit{th}}$$ parameter means main effect plus higher-order effect due to interactions of the $$i^{\mathit{th}}$$ parameter. In this study, the interaction effects of parameters on the output (QoI) are negligible, therefore the main effects are used for factor fixing and ranking of bifurcation parameters.

3.2. Algorithm to compute sensitivity indices

In this section, a detailed working algorithm is presented to compute the main effect, $$S_i$$ using the Monte Carlo simulations, we follow the steps, given in [1, 15].

1. Generate a random numbers matrix of row dimension $$2K$$ and column length $$\mathrm{N}$$ (the sample size) and split into two independent sampling matrices, $$A\left(\mathrm{N},K\right)$$ and $$B\left(\mathrm{N},K\right)$$ by using LHS. Where, $$K$$ is the number of uncertain model parameters.

   $$A\left(\mathrm{N},K\right)=\left[\begin{array}{llll}{x}_{11}& x_{12}& \dots & x_{1K}\\ x_{21}& x_{22}& \dots & x_{2K}\\ \dots & \dots & \dots & \dots \\ x_{\mathrm{N}1}& x_{\mathrm{N}2}& \dots & x_{\mathrm{N}K}\end{array}\right]$$ E13
   $$B\left(\mathrm{N},K\right)=\left[\begin{array}{llll}{x}_{1\left(K+1\right)}& x_{1\left(K+2\right)}& \dots & x_{1\left(2K\right)}\\ x_{2\left(K+1\right)}& x_{2\left(K+2\right)}& \dots & x_{2\left(2K\right)}\\ \dots & \dots & \dots & \dots \\ x_{\mathrm{N}\left(K+1\right)}& x_{\mathrm{N}\left(K+2\right)}& \dots & x_{\mathrm{N}\left(2K\right)}\end{array}\right]$$ E14

2. Define matrix $$C_i$$ , which is matrix A except the $$i^{\mathit{th}}$$ column of matrix B.

   $$C_i\left(\mathrm{N},K\right)=\left[\begin{array}{llllll}{x}_{11}& x_{12}& \dots & x_1\left(K+i\right)& \dots & x_{1K}\\ x_{21}& x_{22}& \dots & x_2\left(K+i\right)& \dots & x_{2K}\\ \dots & \dots & \dots & \dots & \dots & \dots \\ x_{\mathrm{N}1}& x_{\mathrm{N}2}& \dots & x_{\mathrm{N}}\left(K+i\right)& \dots & x_{\mathrm{N}K}\end{array}\right]$$ E15

3. Compute and save model runs for all parameter spaces using matrices $$A$$ , $$B$$ and $$C_i$$ i.e. $$Y_A\left(t,T_s,\mathrm{N}\right)=f\left(A\right),Y_B\left(t,T_s,\mathrm{N}\right)=f\left(B\right)\text{and}{Y}_{C_i}\left(t,T_s,\mathrm{N},K\right)=f\left(C_i\right)$$ , where, $$t$$ are the time points for one heart beat with period $$t_p=0.8s$$ , $$T_s$$ represents the state variables (pressure and flow time series at six locations of arm artery $$\left({\mathrm{N}}_{T_s}=15\right)$$ and $$\mathrm{N}$$ is the total number of model runs $$\left(\mathrm{N}=4000\right)$$ .

4. For the time dependent model outputs, we compute the time dependent main sensitivity index, of each parameter at each time-point of the pressure and flow waves, using the estimator offered by Jansen [15, 16, 17].

   $$S_{t_i}=\frac{V_i}{V}=\frac{V_{xi}\left(E_{x\sim i}\left(Y|x_i\right)\right)}{V}=\frac{V-\frac{1}{2\mathrm{N}}\sum _{n=1}^{\mathrm{N}}{\left(Y_B^{\left(n\right)}-Y_{C_i}^{\left(n\right)}\right)}^2}{V}=1-\frac{1}{2\mathrm{N}\times V}\sum _{n=1}^{\mathrm{N}}{\left(Y_B^{\left(n\right)}-Y_{C_i}^{\left(n\right)}\right)}^2$$ E16

   where,

   $$V=\frac{1}{\mathrm{N}}\sum _{n=1}^{\mathrm{N}}{\left(Y_B^{\left(n\right)}\right)}^2-E^2$$ E17

   and

   $$E={\left(\frac{1}{\mathrm{N}}\sum _{n=1}^{\mathrm{N}}{Y}_B^{\left(n\right)}\right)}^2$$ E18

   The total variance $$\left(V\right)$$ and the expectation $$\left(E\right)$$ are also calculated at each time-point of pressure and flow waves with respect to each parameter.

5. Finally, the main effect, Si of each parameter on the state variables is calculated.

   $$S_i=\frac{1}{{\mathrm{N}}_{T_s}}\frac{1}{{\mathrm{N}}_t}\sum _{j=1}^{{\mathrm{N}}_{T_s}}\sum _{t=0}^{{\mathrm{N}}_t}{S}_{t_i}\left(t,j,t\right),i=1,2,\dots ,K$$ E19

   In Eq. (19), $${\mathrm{N}}_{T_s}$$ is the number of output variables (pressure and flow time series at all locations) and $${\mathrm{N}}_t$$ is the number of time-points [12].

3.3. Input parameters distribution

The results of the UA and SA are greatly affected by the choice of input parameters distributions. In principle, the parameters distributions should be estimated using medical data. Unfortunately, the medical data is not easy to obtained. The input parameters distributions could be chosen according to the expert opinion or using the data from the literature. Due to limited data availability, here in this work the input parameters are randomized within ±10% range of their base (nominal) values using Latin hypercube sampling (LHS).

3.4. Convergence of sensitivity indices

The method of Sobol requires $$\mathrm{N}\left(K+2\right)$$ number of model simulations to compute $$S_i$$ . The main effect, $$S_i$$ is computed for $$\mathrm{N}$$ = [500, 1000, 2000, 3000, 4000] model runs. It is observed that, when the total number of simulations run $$\mathrm{N}$$ increases from 3000 then the sensitivity indices ( $$S_i$$ ) become stable [18]. Therefore, the minimum number of simulations for each parameter to achieve convergence of sensitivity indices is around 3000.

4. Results and discussion

In this section, the sensitivity results based on main effect $$S_i$$ are presented. In order to calculate sensitivity time series, the method of Sobol is applied on each time point of the output QoI i.e. pressure and flow waves at each location of the arm arteries. For each parameter, there are two sensitivity time series at each segment of the arm arteries, one for the pressure and one for the flow. In total, $$K\times {\mathrm{N}}_{T_s}=45\times 33=1485$$ sensitivity time series are obtained. In order to represent the sensitivity results in a compact way, mean absolute values of each pressure and flow sensitivity time series per parameter is taken. In this way, a matrix of dimension $$45\times 33$$ is acquired, where each entry of the matrix represents the mean absolute values of pressure and flow sensitivity time series per parameter, see Figure 3. The numbers in the boxes show the impact (%) on the output (pressure and flow) when input parameters $$\left(R_i,C_i,L_i\right)$$ are randomized within the feasible ranges of ±10%. The parameters having main effect, $$S_i>10\%$$ on output QoI are not shown in the Figure 3. Each row in Figure 3 represents the ranking of influential (bifurcation) parameters. For convenience, the electrical parameters $$\left(R_i,C_i,L_i\right),i=1,2,3,\dots ,15$$ that have impact greater than 10% on pressure and flow are considered as bifurcation parameters which further can be used in bifurcation analysis. For example, for pressure at node-2, $$L_1$$ and $$L_2$$ are the bifurcation parameters, see in Figure 4 (top). Whereas, for flow at node -2, $$L_1$$ and $$L_2$$ are considered as bifurcation parameters, see Figure 4 (bottom).

Figure 3.

Main effect sensitivity of 45-electrical parameters $$\left(R_i,C_i,L_i\right)$$ on pressure and flow time series $$\left(p_i,q_i\right)$$ at 15-segments $$\left(N_1,N_2,\dots ,N_{15}\right)$$ of the arm arteries with number of simulations run per parameter, N is 4000. In total N(K + 2) = 4000(45 + 2) = 188,000 = 0.188 million of simulations run are required to compute the main and total sensitivity indices. The total time taken to compute the sensitivity indices is approximately 3 hours.

Figure 4.

Ranking of bifurcation parameters $$\left(R_i,C_i,L_i\right)$$ in complete arm arteries for pressure (top) and flow (bottom) at node-2. It can be clearly seen that, $$L_1$$ , $$L_2$$ and $$L_1$$ , $$L_3$$ are considered as bifurcation parameters for pressure and flow at node-2 respectively.

In a similar fashion, each row of Figure 3 represents the ranking of bifurcation parameters which further can be used in bifurcation analysis. The parameters which have main effect $$S_i<10\%$$ can be exempted from the bifurcation analysis. The criteria for factor fixing vary from problem to problem.

5. Conclusion

In this chapter, a 5-step procedure of global sensitivity analysis is presented to identify the bifurcation parameters in a lumped-parameter model of the arm arteries. Moreover, the proposed procedure can be applied on any morphology or structure of the systemic circulation (carotid bifurcation, aorta or complete systemic circulation). The results of sensitivity analysis are useful to identify and rank the bifurcation parameters, as well as help which parameters could be exempted from the bifurcation analysis. In this particular example of the arm arteries, 23 out of 45 parameters can be excluded from the bifurcation analysis. Whereas, 22 identified as bifurcation parameters, which further can be used/studied in the bifurcation analysis.