A Priori Identifiability Analysis of Cardiovascular Models.

HHS Public Access Author manuscript Author Manuscript

Cardiovasc Eng Technol. Author manuscript; available in PMC 2015 December 30. Published in final edited form as: Cardiovasc Eng Technol. 2013 December ; 4(4): 500–512. doi:10.1007/s13239-013-0157-3.

A Priori Identifiability Analysis of Cardiovascular Models Jonathan A. Kirk1, Maria P. Saccomani2, and Sanjeev G. Shroff1,3 1Department

of Bioengineering, University of Pittsburgh, Pittsburgh, PA, USA

2Department

of Information Engineering, University of Padova, Padua, Italy

3307

Center for Bioengineering, 300 Technology Drive, Pittsburgh, PA 15219, USA

Author Manuscript

Abstract

Author Manuscript

Model parameters, estimated from experimentally measured data, can provide insight into biological processes that are not experimentally measurable. Whether this optimized parameter set is a physiologically relevant complement to the experimentally measured data, however, depends on the optimized parameter set being unique, a model property known as a priori global identifiability. However, a priori identifiability analysis is not common practice in the biological world, due to the lack of easy-to-use tools. Here we present a program, Differential Algebra for Identifiability of Systems (DAISY), that facilitates identifiability analysis. We applied DAISY to several cardiovascular models: systemic arterial circulation (Windkessel, T-Tube) and cardiac muscle contraction (complex stiffness, crossbridge cycling-based). All models were globally identifiable except the T-Tube model. In this instance, DAISY was able to provide insight into making the model identifiable. We applied numerical parameter optimization techniques to estimate unknown parameters in a model DAISY found globally identifiable. While all the parameters could be accurately estimated, a sensitivity analysis was first necessary to identify the required experimental data. Global identifiability is a prerequisite for numerical parameter optimization, and in a variety of cardiovascular models, DAISY provided a reliable, fast, and simple platform to provide this identifiability analysis.

Keywords Identifiability analysis; Systemic arterial circulation models; Cardiac muscle contraction models; Numerical parameter estimation

Author Manuscript

INTRODUCTION With the availability of affordable and powerful computers and high-level programming languages such as MatLab, the complexity of mathematical models of biological systems has grown significantly. Model parameters, estimated from experimentally measured data (selected inputs and outputs), are often used to provide mechanistic insights into the

Address correspondence to Sanjeev G. Shroff, 307 Center for Bioengineering, 300 Technology Drive, Pittsburgh, PA 15219, USA. [email protected]. ELECTRONIC SUPPLEMENTARY MATERIAL The online version of this article (doi:10.1007/s13239-013-0157-3) contains supplementary material, which is available to authorized users.

Kirk et al.

Page 2

Author Manuscript Author Manuscript

biological processes that are not directly measurable. Specifically, model parameters are iteratively adjusted to minimize the difference between model-based and experimentally measured outputs and some statistical criterion is used to judge the optimality of parameter estimates. It is assumed that this optimized parameter set is a significant and physiologically relevant complement to the experimentally measured inputs and outputs. The reliability of this assumption, however, often depends on whether the parameter set is unique. The a priori model identifiability analysis can theoretically address this question. Specifically, this analysis checks the uniqueness of model parameters based on the model structure (equation system) and ideal inputs and outputs (i.e., noise-free and no limitations on the information content). There are three possible results from this analysis: (1) globally identifiable (i.e., it has a unique set of model parameters), (2) locally identifiable (multiple, but finite, parameter sets), or (3) non-identifiable (infinite number of parameter sets). In the case of a nonidentifiable set of model parameters, the analysis can be helpful in providing guidance for modifying the experiments or model structure in order to reach global identifiability.

Author Manuscript

The a priori identifiability analysis should be implemented before conducting experiments and performing a model-based analysis of experimentally measured input–output data. However, this is not a common practice in the biological world, primarily due to the lack of easy-to-use analytical tools for performing this type of analysis. For non-linear models, no standard algorithm exists for testing a priori global identifiability. Several approaches to identifiability have been developed (see Bellu et al.3 for examples) to analyze linear models, but to the uninitiated, they can be very difficult to apply. With the aid of computer algebra packages, i.e. Maple, Mathematica, it is possible to implement symbolic mathematical calculations required for a variety of identifiability analysis approaches. However, in general, the required mathematical background may not be available to investigators who use experimental data to identify parameters in specific dynamic models, such as cardiologists or physiologists. Here we present a software tool which can be used without requiring high-level programming languages, mathematics, or computer algebra as prerequisites. A recently developed program that uses a differential algebra method, named Differential Algebra for Identifiability of Systems (DAISY), greatly facilitates this process for ordinary linear and/or nonlinear differential equation-based models.3,16,21 It has been shown that the differential algebra algorithm implemented by DAISY improves the efficiency of the previous differential algebra algorithms models and enlarges their applicability domain.24

Author Manuscript

DAISY can be directly applied to a wide-variety of models incorporating polynomial or rational equations. More generally, it can be also applied to test identifiability of all models described by expressions that can be represented as solutions of differential equations with a rational right hand side by introducing new state variables,1 including exponential functions with integer exponent and models characterized by time-varying parameters. However, there are some functions (PDEs, time delays, non-integer exponentials) which DAISY cannot analyze in an automated way. In some of these cases, however, approximations can be utilized.

Cardiovasc Eng Technol. Author manuscript; available in PMC 2015 December 30.

Kirk et al.

Page 3

Author Manuscript

The primary goal of the present study was to examine the utility of a DAISY-based a priori identifiability analysis in the context of several commonly used cardiovascular models. Specifically, we wish to address the following two questions: (1) Can DAISY provide information regarding a priori identifiability of these models, given specific model structures and experimental data (i.e., inputs and outputs)? (2) Can DAISY provide insights into what makes a given model non-identifiable so that appropriate actions can be taken to make it identifiable (e.g., fixing certain model parameters, additional inputs/outputs, change in model structure, etc.)?

Author Manuscript

The globally identifiable outcome of the a priori identifiability analysis may not be sufficient to guarantee a reliable estimation of the model parameters from real (experimental) input/output data wherein, unlike ideal data, noise is present and information content is limited. Therefore, the secondary goal of the present study was to examine this issue of model parameter estimation using real data for cardiovascular models that were judged to be globally identifiable.

METHODS A Priori Model Identifiability We assume the model is defined by a system of differential equations: (1)

where x is an n-dimensional state variable vector, p is an i-dimensional parameter vector, and u is an m-dimensional input vector. We also assume there is at least one output equation:

Author Manuscript

(2)

where y is a r-dimensional output vector. Both f and g are assumed to be polynomial or rational functions. The identifiability analysis is in the context of a specified model structure and experimental data (i.e., specified inputs and outputs). In addition, it assumes ideal inputs and outputs (i.e., noise-free and without any constraints on the information content). The model is globally (uniquely) identifiable if there exists a unique parameter set, p = P, satisfying Eqs. (1)–(2).

Author Manuscript

Methods have been proposed to check identifiability of linear and nonlinear dynamical models.13,16,21,29 Differential algebra techniques can be used to study identifiability in a certain class of dynamical systems, i.e. systems described by polynomial or rational differential equations.17 A brief description of the differential algebra method recently implemented by the software DAISY is presented in this section. For further details the reader is referred to.3,24 Consider the set of equations (Eqs. (1)–(2)) as being defined by a system of n + r (number of state variables plus outputs) differential polynomials with rational coefficients and with variables x, y and u (state variable, output, and input, respectively). The goal is to calculate


Kirk et al.

Page 4

Author Manuscript

the characteristic set of the model, which is the “minimal” set of differential polynomials that, when set to zero, has the same solutions as the original model (Eqs. (1)–(2)). To accomplish this, one has to use the Ritt’s algorithm, which is analogous to the Gauss elimination algorithm used to solve linear algebraic equations. This requires the introduction of a ranking among the variables; in particular, the higher ranked variables are eliminated first. In the context of parameter identifiability, the unknown state variables and their derivatives are ranked highest so they are preferentially eliminated. Then the polynomials defining the system (Eqs. (1)–(2)) are reduced using Ritt’s pseudo-division algorithm. Once the system can no longer be reduced, the characteristic set has been obtained.

Author Manuscript

It has been shown that if the original dynamic system is in its state-space form, as in Eq. (1), the characteristic set has a diagonal structure, and the first r (number of outputs) polynomials of the characteristic set contain only input and output variables.21,24 Therefore, they are obtained after elimination of the state variables x from the system (Eqs. (1)–(2)). The corresponding system of differential polynomial equations: (3)

Author Manuscript

where the number of derivatives are fixed from theoretical considerations on the system order, thus represents the so-called input–output relation of the model. This describes all input–output pairs which satisfy the system equations (Eqs. (1)–(2)). These r equations are the basic tool used in the identifiability analysis. Equation (3) is therefore a function of only the model inputs, outputs, and their derivatives, which, under the a priori identifiability hypotheses, are all theoretically known variables. Hence its coefficients in the parameter vector, p, are known. After a suitable normalization to make all the polynomials monic in order to uniquely fix the coefficients, these are extracted to form the exhaustive summary of the model. This is a nonlinear algebraic system in the unknown parameters (p) of the model. The Buchberger algorithm is applied to solve it. In particular, this computer algorithm calculates the Gröbner basis of the system which provides the number of solutions for each unknown parameter.18 If all the parameters have only one solution, then the model is a priori globally identifiable. If there are a finite number of possible parameter sets then the model is locally identifiable, and if there are an infinite number of parameter sets the model is non-identifiable.

Author Manuscript

The theoretical details of this analysis are explained in detail elsewhere.1,3 The major steps of the differential algebra algorithm used to determine whether a model is globally identifiable are easiest to describe through a simple example. A Priori Identifiability Analysis: An Example Consider a two compartment model (see Fig. 1), with an input u(t) into compartment 1 (x1), a transition rate constant k1 from compartment 1 to compartment 2 (x2) and a rate constant k2 representing a sink from compartment 2. The model output, y, is the measurement in


Kirk et al.

Page 5

Author Manuscript

compartment 2, of x2. The model is described by the following two differential equations and one output equation. (4)

(5)

(6)

Author Manuscript

In the following, we give an idea how the pseudo-division algorithm works, without stating all the formal definitions and theorems on which it is based. We will begin by assigning the standard ranking to all of the variables. The state variables will be ranked the highest and thus eliminated first, because the goal is to reduce the equations to an input–output relationship. (7)

The highest ranked variable in each of the polynomials defining the model (Eqs. (4)–(6)) is identified (called the “leader”). Then the equations are arranged according to the highest ranked leader. This yields: (8)

Author Manuscript

(9)

(10)

The first step is to reduce the highest ordered polynomial, Eq. (10), by using Eq. (8) and its derivative. This step reduces the ranking of the leader, and transforms Eq. (10) into: (11)

The system is now:

Author Manuscript

(12)

(13)

(14)


Kirk et al.

Page 6

Author Manuscript

The same technique is used again, where Equation (14) now contains the highest ranked variable, using Eq. (12) and its derivative are used to replace . When this is done, Eq. (14) becomes: (15)

The system is now: (16)

(17)

Author Manuscript

(18)

This system cannot be further reduced, because the highest ranked variables, x1 and x2, cannot be eliminated from Eqs. (17) and (18), respectively. Therefore, Eqs. (16)–(18) represent the model’s characteristic set. In this case the number of outputs is 1 thus the first polynomial (Eq. (16)) represents the input–output relation of the model. This is already monic, thus its coefficients can be extracted to form the exhaustive summary: (19)

Author Manuscript

The next step is to determine if k1 and k2 can be uniquely identified from the exhaustive summary. By equating the polynomials (Eq. (19)) to known symbolic values α, β and γ, we obtain a nonlinear equation system in the unknown k1 and k2 and the Buchberger algorithm is applied to obtain the Gröbner basis. For this simple example, it is possible to do this analysis by hand: (20)

(21)

(22)

Author Manuscript

Clearly, this system can be solved for k1 = γ and k2 = β/γ. This technique has been automated in DAISY. Therefore, for this example, DAISY would report that the model described by Eqs. (4)–(6) is globally identifiable. The calculation of the Gröbner basis, with symbolic coefficients, can become computationally intensive with even minor increases in model complexity, making it difficult or impossible to solve. In order to assuage this problem, DAISY uses a pseudorandomly selected numerical value instead of the symbolic one, to simplify the symbolic


Kirk et al.

Page 7

Author Manuscript

algebra, and to make it possible to perform identifiability analysis for complex models; this is theoretically justified in Bellu et al.3 In the case of a locally identifiable, or non-identifiable model, initial conditions can be provided if they are known.24 In some cases, supplying known initial conditions can make the model globally identifiable. Moreover, in a globally identifiable model, it is possible to test the identifiability of a set of unknown initial conditions. DAISY (version 1.5) was written by us (M. P. Saccomani), and others (G. Bellu,, S. Audoly and L. D’Angio), and is copyrighted by University of Cagliari and University of Padova, Italy. The program runs in Reduce (version 3.8, Konrad-Zuse-Zentrum für Informationstechnik Berlin (ZIB), Berlin, Germany), an algebraic computer program which is available for free (http://reduce-algebra.sourceforge.net/).

Author Manuscript

Cardiovascular Models The cardiovascular system is a common target for mathematical modeling and consists of two subsystems: the heart acting as a pump, and an organization of vessels acting as conduits. There are two classes of models commonly used in this area: systemic arterial circulation models, and cardiac muscle contraction models. Both types of models can predict or explain the behavior of the cardiovascular system.

Author Manuscript

Oftentimes, the biological processes of interest cannot be directly measured. For example, there currently exists no easy method for directly measuring the binding constant for calcium to troponin C under in vivo conditions, a key step in myofilament activation, or for determining arterial compliance in a live animal. Model parameters, estimated using experimentally measured data, are often used to provide mechanistic insights into the biological processes that are not directly measurable. Below, several cardiovascular models are presented that will be analyzed for a priori identifiability. Systemic Arterial Circulation—Windkessel Models Aortic pressure and flow are linked via the mechanical properties of the systemic arterial circulation. Lumped parameter models with analogous electrical circuit representation are extensively used to characterize systemic arterial circulation. One of the simplest forms is the three-element Windkessel model28,30; consisting of peripheral arterial resistance (RP), total arterial compliance (C) and characteristic impedance (r), as illustrated in Fig. 2a. The model is described by one differential equation, Eq. (23), and one output equation, Eq. (24).

Author Manuscript

(23)

(24)

where PC is a state variable for pressure (voltage) across the capacitor, Qao is the aortic flow (model input) and Pao is the aortic pressure (model output). By following the procedure reported in the previous section, the state variable PC and its derivative are ranked the


Kirk et al.

Page 8

Author Manuscript

highest in the DAISY analysis:

, thus the leader

of Eq. (23) is , and the leader of Eq. (24) is PC. Arterial total vascular resistance (TVR) can be directly estimated from the input and output (TVR = mean Pao/mean Qao = RP + r); so there are only two free parameters to be estimated: r and C (or Rp and C). While the three-element Windkessel model captures many of the major characteristics of the aortic pressure and flow relationship, greater fidelity comes from including additional elements. For example, the five-element Windkessel28 includes the effects of blood inertance,27 and has a more complicated structure (Fig. 2b). The model is described by three differential equations and one output equation. (25)

Author Manuscript

(26)

(27)

(28)

Author Manuscript

where PC1, QL, and PC2 are the state variables, Pao is the output, and Qao is the input. There are four free parameters to be estimated: r, C1, C2, and L. RP is equal to TVR in this model; so this resistance does not have to be estimated. Systemic Arterial Circulation—T-Tube Model Windkessel models have lumped parameters that do not account for the spatially distributed nature of the systemic arterial circulation and consequently, phenomena such as wave reflections. The T-Tube model consists of two asymmetrical wave transmission tubes, one representing blood flow to the head and one representing blood flow to the body, both of which originate at the heart,5,6,9 as shown in Fig. 2c. At the end of each tube there is a complex load, represented by a three-element Windkessel model. There are two linked differential equations, representing the pressure in the body portion, Pcb, and in the head, Pch.

Author Manuscript

(29)

(30)


Kirk et al.

Page 9

Author Manuscript

(31)

Pao and Qao are the model input and output, respectively, and . There are a total of ten parameters: Zcb, Zch, τb, τh, Clb, Clh, Rob, Roh, Rpb, and Rph. This parameters space can be reduced to seven free parameters by enforcing the following three constraints: (1)

, (2)

(high-frequency

Author Manuscript

matching of the body-ward with its terminal load),26 (3) (high-frequency 26 matching of the head-ward tube with its terminal load). Thus, the seven free parameters are: τb, τh, Clb, Clh, Rob, Roh, and Rpb (or Rph). The τ parameters are time delays, which DAISY cannot handle, and so were replaced with second order Taylor series expansions. Shroff et al.26 measured a second output, the blood flow to the body, Qb(t): (32)

where τd is a known time constant (delay between upstroke of Pao and the onset of Qb). The second model output makes it possible to calculate Rpb and Rph, since pressure and flow are known for each load (Rpb = mean Pao/mean Qb and Rph = mean Pao/(mean Qao − mean Qb)), leaving six free parameters to estimate: τb, τh, Clb, Clh, Rob, and Roh.

Author Manuscript

Cardiac Muscle Contraction—Stiffness Model Cardiac muscle responds to changes in length with alterations in active force development; this force–length relationship (FLR) represents a dynamic process. It is common to alter muscle length sinusoidally, and observe the resultant sinusoidal force response. The ratio of the force response and length perturbation at each sinusoidal frequency represents the complex stiffness (i.e., having magnitude and phase) at that particular frequency. The complex stiffness over a wide range of frequencies (experimentally obtained by a length “chirp”, a sinusoid of continuously increasing frequency) is a useful indicator of a muscle’s dynamic FLR.

Author Manuscript

The stiffness model was developed to analyze chirp data from skinned fibers in the time domain, and is based on system of two force producing dynamic processes: recruitment and distortion of crossbridges. The model has been rigorously validated,10 and used for parameter estimation,11,12 but has never been tested for identifiability. The model is described by two differential equations: Eqs. (33) and (34) represent the recruitment and distortion dynamics, respectively, and Eq. (35) is the output equation: (33)


Kirk et al.

Page 10

Author Manuscript

(34)

(35)

where η(t) is the recruitment variable, x(t) is the distortion variable, ΔL is the length change (model input), and ΔF is the force change (model output). There are four free parameters to be estimated: E0, E∞, b, and c. Cardiac Muscle Contraction—Three-State Model

Author Manuscript

The sarcomere is comprised of two main units, the thick filament and the thin filament.20 The thick filament contains the myosin heads, which can bind to actin on the thin filament to form a crossbridge, the force-generating structure. Tropomyosin blocks the interaction of myosin and actin until calcium binds to troponin, which causes a conformational shift in tropomyosin which allows the formation of a cross-bridge. Calcium activates the sarcomere to produce force, which drives the ventricle to generate pressure. There are few experimental techniques to reliably measure myofilament properties, so their measurement is handled indirectly, through identifying model parameter values. There are several models of myofilament function that can be used for this approach, often taking the form of compartment or state models. The first of these examined here was developed to analyze the effects of acidosis,4 and is shown in Fig. 3a. It is described by the following three differential equations,

Author Manuscript

(36)

(37)

(38)

(39)

Author Manuscript

On the right hand side of the equations, d is a parameter, not a derivative. The three states: Xco (cycling unattached crossbridges), Xcn (cycling, attached, non-force-producing crossbridges), and Xcf (cycling, attached, force-producing crossbridges), along with Xnc (non-cycling crossbridges) sum up to XT (total number of crossbridges, which is a constant parameter, XT = Xco + Xcn + Xcf + Xnc). The input is the intracellular calcium concentration, [Ca] (either a time-dependent transient, or a steady state value) and the output is F. The parameter S is a scaling parameter that also needs to be estimated. Thus, there are eight free parameters to estimate: koff, kon, f, d, g, h, XT and S.


Kirk et al.

Page 11

Author Manuscript

Initial conditions for the state variables were provided to DAISY as well. These were provided for the three-state model because they define a structural element by restricting the total number of cycling cross bridges. The model can also include a feedback algebraic equation, which will decrease koff when the value of Xcf, force producing crossbridges, gets larger. This feedback attempts to describe the concept of cooperativity in the model. The feedback equation is: (40)

Author Manuscript

where koff,0 is the baseline value of koff, and β is a cooperatively parameter. This removes one free parameter (koff), but introduces two, increasing the number of free parameters to nine. Cardiac Muscle Contraction—Four-State Model The four-state model is a more complicated compartment model, and has been utilized by us and others.2,8,19,22,25 A schematic of the four-state model is shown in Fig. 3b. It is defined by three differential equations for each of the species ([A], actin; [M], myosin; [Ca–A], calcium bound to troponin/actin; [Ca–A–M], force generating crossbridge in the presence of bound calcium; [A–M], force generating crossbridge without bound calcium) and two feedback equations (force generating states feedback on time-varying parameters k1 and f), and one output equation: (41)

Author Manuscript

(42)

(43)

(44)

Author Manuscript

(45)

(46)


Kirk et al.

Page 12

Author Manuscript

(47)

(48)

The model input is [Ca](t), and the model output is F(t). The exponents on the feedback equations (Eqs. (46) and (47)), γ1 and γf, are 0.5 and 2, respectively. There are ten free parameters: α1, β1, k2, k3, k4, αf, βf, d, g and S. Identifiability analysis has also been performed without the feedback equations, in which case there are eight free parameters (removing α1, β1, αf, and βf but introducing the parameters k1 and f). As with the three-state model, the total number of cycling crossbridges was restricted by providing initial conditions.

Author Manuscript

DAISY cannot handle exponents less than 1, so there is an issue with Eq. (46) when γ1 is equal to 0.5. In order to convert the equations in polynomial form, we follow the strategy of introducing a new state variable, see Appendix for the details.

RESULTS AND DISCUSSION

Author Manuscript

Using DAISY, identifiability was determined for six cardiovascular models. The DAISY output includes: a reiteration of the model (equations, inputs, outputs, variables), the variable ranking (see Eq. (7)), the characteristic set, the randomly selected parameter vector, the exhaustive summary, the parameter solutions, and a summary of model identifiability (Globally-, Locally-, or Non-Identifiable). An example DAISY input and output file is shown in the Appendix. The exhaustive summary is often too complex to interpret (the example is relatively simple), but the parameter solutions can yield additional information. The results obtained are summarized in Table 1. Globally Identifiable Models The three-element and five-element Windkessel models are relatively simple, and DAISY finds these models to be a priori globally identifiable. These models, and a variety of variations, have been used extensively. They have been fit to experimental data using standard optimization techniques for decades, and provide excellent reproduction of the data, as well as physiologically relevant parameter estimates.7 While a good fit does not necessarily imply an identifiable model, the extensive use of these models highly suggest they are globally identifiable, and so it is encouraging that DAISY finds them to be so.

Author Manuscript

The dynamic muscle stiffness model was also found to be globally identifiable. In addition, both the three-state and four-state models were globally identifiable. This was true regardless of the presence or form of the feedback equations. One of the keys to an identifiable state model is distinctiveness of the pathways described by each parameter, as was true of both of the state models presented here. Consider a state model where one of the states had two sinks, and the rate of loss through each sink is controlled by a separate parameter, but the only measurement is of total loss (sum of the two sinks). There is no information in the model that allows it to distinguish between these two sinks (and by


Kirk et al.

Page 13

Author Manuscript

extension, the associated parameters); therefore this model would be non-identifiable. Thus, an identifiable state-model should have distinct pathways, and avoid simple parallel redundancy. Again, the three-state and four-state models possessed both of these qualities. While most biological systems have some level of redundancy, these often invoke other states and pathways, and thus are distinct and not parallel. In some cases, the model must be simplified. By combining, simplifying, or bounding parameters, it may be possible to turn a non-identifiable model into an identifiable one. Non-Identifiable Model

Author Manuscript

The T-Tube model with one input (Pao(t)), one output (Qao(t)), and seven unknown parameters was non-identifiable. Daisy provides the uniquely identifiable parameterization, i.e. the aggregates of model parameters that are globally identifiable. If these parameters are of sole interest to the investigator, they can proceed with their numerical estimation from the experimental data. As mentioned previously, DAISY can be used to determine what additional information can be gathered in order to make the model identifiable. In particular, DAISY output file reported the following Gröbner basis for a randomly preselected parameter values: Rpb = 21, Rob = 18, Roh = 27, Clb = 25, Clh = 23, τb = 8, τh = 10 (numerical values were used, instead of symbolic parameters, to speed up the computation process): Rpb = 21, Rob = 18, Roh = 27, Clb = (5/2) * τh, Clh = (23/10) * τh, τb = (4/5) * τh.

Author Manuscript

It can be seen that while parameters Rpb, Roh and Rob were uniquely identifiable, the other four parameters, Clb, Clh, τb, and τh were not. τh has an infinite number of solutions, and Clb, Clh, and τb are all linearly related to τh. Thus, all four of these parameters are related, and if one of these parameters were assumed to be known, the model would be globally identifiable. The non-uniqueness in the estimation of these parameters indicates that one cannot uniquely differentiate between the properties of the body and head portions of the model and suggests that some extra information might make the model globally identifiable. If one or more of these parameters could be identified via direct experimental measure, this would be ideal. The additional information can come from a variety of other sources, such as fixing a parameter value based on a literature search or a logical assumption or having additional experimental data.

Author Manuscript

Shroff et al.26 used the T-tube model to fit experimental data with one input and one output, but the parameter estimation was poor. It was possible, in their experimental system, to measure a second output, the flow to the body portion of the model, Qb(t). The addition of this output variable into the model reduces the number of unknown parameters by one, since both Rph and Rpb can now be calculated, bringing the number of unknown parameters to six. In addition, this extra measurement made it possible to uniquely differentiate between the body-end and head-end components of the model. Indeed, with the addition of the second output measurement, DAISY-based analysis indicated that the T-Tube model was globally identifiable. This underscores the importance and versatility of a priori identifiability analysis when designing experiments.


Kirk et al.

Page 14

Numerical Parameter Estimation


Once the a priori identifiability analysis is complete, however, the process is not over, as there are other issues that can arise. DAISY’s identifiability analysis determines whether there is a unique parameter set satisfying the given model-experiment structure,14 not whether that parameter set is practically identifiable from a typical set of experimental data. Therefore, we will continue with numerical parameter estimation for a model which has been verified by DAISY to be globally identifiable, specifically, the four-state model of cardiac muscle contraction. Parameter estimation was performed using the Matlab function lsqcurvefit (which uses the Levenberg–Marquadt and Newton methods). Each parameter was bounded by zero as a minimum, but given no maximum bound. The initial conditions represent the concentrations (per unit volume of sarcoplasm) of each state prior to introduction of free calcium2,23: [A]t=0 = 70 μM, [M]t=0 = 20 μM, [Ca-A-M]t=0 = 0 μM (these known initial conditions were supplied to DAISY, but because the model was found to be globally identifiable, they were unnecessary). The four-state model can be used to fit two types of experimental data: dynamic data (timevarying force and calcium transients from twitching muscle or myocytes), or steady-state data (force-calcium relationships from membrane permeabalized preparations that are exposed to constant calcium levels). To gauge how well the model parameters can be estimated, their values must be known beforehand. Therefore, a synthetic data set was utilized, generated using the four-state model (Eqs. (41)–(48)) with an assumed parameter set and synthetic input calcium. The dynamic synthetic input (calcium transient) was generated using:

Author Manuscript

(49)

where A = 1.739 μM, B = 205.89 s−1, and C = 30.52 s−1. For the steady-state synthetic input, constant levels of calcium were used as input. The parameter set used to generate the synthetic data (both dynamic and steady-state) is listed in Table 2. The synthetic timevarying force waveform is shown in Fig. 4 (solid line).

Author Manuscript

Numerical optimization was used to determine whether it was possible to estimate the unknown parameter set from this synthetic data. The second column of Table 2 shows the parameter values estimated from the dynamic synthetic data alone (and not steady-state). The ± values in Table 2 are the 95% confidence limits of estimation calculated from the Jacobian matrix (not to confuse them with SE calculated from experimental observations). While, some of the parameters were accurately estimated, and associated with reasonable 95% confidence limit values (Table 2, see values highlighted in gray), the majority of the parameters had large confidence limits, indicating they could not be accurately estimated. Therefore, despite a good overall fit (high R2, Table 2, Fig. 4); the actual parameters could not be uniquely identified. This represents an example of a case where there is a good fit, but poor parameter estimation – precisely the reason DAISY is important. However, we know that this model can produce unique parameter sets, as it was globally identifiable. A common practice to improve parameters estimation is to use more data.


Kirk et al.

Page 15

Author Manuscript

While the model is normally fit to dynamic data alone, numerical parameter optimization using synthetic dynamic and steady-state data simultaneously resulted in accurate estimation of all model parameters (Table 2). The standard error values were several orders of magnitude smaller than those associated with parameters fit to only the dynamic data. As stated earlier, DAISY assumes that the input/output pair has no restraints on informational content. It is possible that some of the parameters were insensitive to changes in either the dynamic or steady state data. To determine if this was the case, we performed a

Author Manuscript

sensitivity analysis on the model. Sensitivity was calculated as , where is the sensitivity of Y to X, and is dimensionless. Table 3 shows the effect of changes in each parameter on several calculated indices from the dynamic and steady-state force data (model output). The mean and SEM sensitivity values in Table 3 were calculated from the sensitivity values at 50, 75, 125, and 150% of the baseline value of the parameter.

Author Manuscript

The indices from the dynamic data (Fdev, Trise, Trelax) were most sensitive to changes in k3, g, and d (see Table 3), and these parameters were associated with reasonable standard error values when fit to the synthetic dynamic data alone (Table 2, see values highlighted in gray). These indices were less sensitive to changes in the other parameters, which had comparatively larger standard errors. The exception to this is β1 and βf, which have a large effect on Fdev but are not well estimated from the synthetic dynamic data alone. Looking at Table 3, the effect of changes in β1 and βf on all three indices of the dynamic data is essentially identical (large effect on βdev with minimal effects on Trise and Trelax). This suggests that while the dynamic (time-varying) data are sensitive to changes in these two parameters, their effects cannot be separated from each other, making it difficult to estimate them. However, β1 and βf are differentiated in the indices from the steady-state data (Fmax, pCa50, nH). Here, β1 has a large effect on nH while βf has a large effect on Fmax (whereas β1 has essentially zero effect on Fmax). Thus, with the addition of the steady-state data, it would be possible to accurately estimate both of these parameters. The estimated parameters were all within 3% of the initial parameter values. This small variation could arise from the error tolerances of the curve fitting algorithm, as the program terminated before the parameter values could be precisely recovered. This result indicates that the input/output pair for the dynamic contraction lacked the necessary informational content to estimate all parameters. To compensate for this, when using the four-state model, it should be fit to dynamic and steady-state data simultaneously to achieve the most accurate estimation of all parameters.

Author Manuscript

This confirms that a positive result from the a priori global identifiability analysis constitutes a necessary, but not sufficient, condition for successful parameter identification from real input–output data. However, even if a model is locally identifiable or even nonidentifiable, forward analysis—in silico experiments—can still be useful. In fact, the identifiability analysis performed by DAISY always provides the globally identifiable parameterization of the model, i.e. the groups of parameters uniquely identifiable, and these can thus be numerically estimated.


Kirk et al.

Page 16

DAISY Limitations

Author Manuscript

While DAISY can be directly applied to a wide-variety of models incorporating polynomial or rational equations and expressions that can be represented as solutions of differential equations with a rational right hand side, there are a variety of functions and models it cannot analyze. Obviously, even if DAISY can be applied to test identifiability of many biological models, its domain of applicability is limited. Models including non-integer exponential functions cannot be analyzed by DAISY. Also, DAISY cannot analyze models described by equations including time delays, in which case some approximations are required. An example is provided here by the analysis of the T-Tube model, where it is shown that approximations can be used effectively.

Author Manuscript

The strategy for describing some non-polynomial systems, exponential or logarithmic functions for example, as solutions of differential equations with a rational right hand side has been proposed.15 This simply requires the introduction of additional state variables, and may be easily integrated into the program in the future. This would improve the breadth of models DAISY can analyze. Even with these limitations, DAISY can be used on a widevariety of complex and useful models.

Author Manuscript

DAISY assumes a noise-free data set. This simply indicates that it is impossible for the program to predict the effect experimental noise might have upon the model. While a small amount of random noise would most likely have little to no effect on the ability to optimize parameters, other sources or a larger magnitude of noise could interfere. For example, if there was a large source of high-frequency noise, it might make it impossible to uniquely identify a parameter that describes very fast events. Another possibility is noise which affects only one measurement, skewing the data for one output measurement but not another. An analysis of the magnitude and source of noise can be helpful. However, this is not an issue where DAISY can provide insight. Lastly, while DAISY is generally very fast, even for complex high order non-linear models, some analyses may fail to terminate successfully. The most computationally demanding step is represented by the Buchberger algorithm (embedded REDUCE function) for solving the system of algebraic nonlinear equations provided by the exhaustive summary. Thus, if the model contains multiple nonlinearities, the algorithm may not successfully terminate due to a lack of memory of the system running the application. Summary

Author Manuscript

We used DAISY, an automated a priori identifiability analysis tool, in the context of several cardiovascular models: systemic arterial circulation (Windkessel, T-Tube) and cardiac muscle contraction models (complex stiffness, crossbridge cycling-based). We found them all to be globally identifiable, except for the T-Tube model with one input and one output. DAISY’s output can provide insight into why the model may be non-identifiable, and can help guide the choice of additional experiments, or suggest that a change in the fundamental model structure is necessary. Global identifiability is not sufficient to guarantee that a model will perform well with actual experimental data, since identifiability assumes noise-free data with no constraints on its information content. We applied numerical optimization


Kirk et al.

Page 17

Author Manuscript

techniques to estimate unknown parameters using typical (synthetic) experimental data and a model which DAISY found to be globally identifiable. We were unable to estimate with good precision all of the unknown model parameters. By using a sensitivity analysis, however, it was found that providing additional data (without changing the model structure) could yield good estimates of model parameters. Identifiability is a prerequisite for numerical parameter optimization, and in a variety of arterial circulation and cardiac muscle models, DAISY provided a reliable, fast, and simple platform to provide this identifiability analysis.

Supplementary Material Refer to Web version on PubMed Central for supplementary material.

Author Manuscript

ACKNOWLEDGMENTS This work was supported by the National Institutes of Health grant T32-HL76124 (S.G.S.) and the McGinnis Endowed Chair Funds (S.G.S.).

REFERENCES


1. Audoly S, Bellu G, D’Angio L, Saccomani MP, Cobelli C. Global identifiability of nonlinear models of biological systems. IEEE Trans. Biomed. Eng. 2001; 48:55–65. [PubMed: 11235592] 2. Baran D, Ogino K, Stennett R, Schnellbacher M, Zwas D, Morgan JP, Burkhoff D. Interrelating of ventricular pressure and intracellular calcium in intact hearts. Am. J. Physiol. 1997; 273:H1509– H1522. [PubMed: 9321844] 3. Bellu G, Saccomani MP, Audoly S, D’Angio L. DAISY: a new software tool to test global identifiability of biological and physiological systems. Comput. Methods Programs Biomed. 2007; 88:52–61. [PubMed: 17707944] 4. Berger DS, Fellner SK, Robinson KA, Vlasica K, Godoy IE, Shroff SG. Disparate effects of three types of extracellular acidosis on left ventricular function. Am. J. Physiol. 1999; 276:H582–H594. [PubMed: 9950860] 5. Burattini R, Campbell KB. Modified asymmetric T-tube model to infer arterial wave reflection at the aortic root. IEEE Trans. Biomed. Eng. 1989; 36:805–814. [PubMed: 2759639] 6. Burattini R, Campbell KB. Effective distributed compliance of the canine descending aorta estimated by modified T-tube model. Am. J. Physiol. 1993; 264:H1977–H1987. [PubMed: 8322928] 7. Burattini R, Di Salvia PO. Development of systemic arterial mechanical properties from infancy to adulthood interpreted by four-element windkessel models. J. Appl. Physiol. 2007; 103:66–79. [PubMed: 17303709] 8. Burkhoff D. Explaining load dependence of ventricular contractile properties with a model of excitation-contraction coupling. J. Mol. Cell. Cardiol. 1994; 26:959–978. [PubMed: 7799451] 9. Campbell KB, Burattini R, Bell DL, Kirkpatrick RD, Knowlen GG. Time-domain formulation of asymmetric T-tube model of arterial system. Am. J. Physiol. 1990; 258:H1761–H1774. [PubMed: 2360669] 10. Campbell KB, Chandra M, Kirkpatrick RD, Slinker BK, Hunter WC. Interpreting cardiac muscle force-length dynamics using a novel functional model. Am. J. Physiol. Heart Circ. Physiol. 2004; 286:H1535–H1545. [PubMed: 15020307] 11. Campbell KB, Wu Y, Simpson AM, Kirkpatrick RD, Shroff SG, Granzier HL, Slinker BK. Dynamic myocardial contractile parameters from left ventricular pressure-volume measurements. Am. J. Physiol. Heart Circ. Physiol. 2005; 289:H114–H130. [PubMed: 15961371]


Kirk et al.

Page 18

Author Manuscript Author Manuscript Author Manuscript Author Manuscript

12. Chandra M, Tschirgi ML, Rajapakse I, Campbell KB. Troponin T modulates sarcomere lengthdependent recruitment of cross-bridges in cardiac muscle. Biophys. J. 2006; 90:2867–2876. [PubMed: 16443664] 13. Chapman MJ, Godfrey KR, Chappell MJ, Evans ND. Structural identifiability for a class of nonlinear compartmental systems using linear/non-linear splitting and symbolic computation. Math. Biosci. 2003; 183:1–14. [PubMed: 12604132] 14. Cobelli C, DiStefano JJ III. Parameter and structural identifiability concepts and ambiguities: a critical review and analysis. Am. J. Physiol. 1980; 239:R7–R24. [PubMed: 7396041] 15. D’Angiò, L.; Saccomani, M.; Audoly, S.; Bellu, G. Identifiability of nonaccessible nonlinear systems. In: Bru, R.; Romero-Vivó, S., editors. Positive Systems. Springer; Berlin/Heidelberg: 2009. p. 269-277. 16. Denis-Vidal L, Joly-Blanchard G. Equivalence and identifiability analysis of uncontrolled nonlinear dynamical systems. Automatica. 2004; 40:287–292. 17. Forsman, K. Linköping Studies in Science and Technology, Department of Electrical Engineering. Linköping University; Linköping, Sweden: 1992. Constructive Commutative Algebra in Nonlinear Control Theory, Dissertation No. 261. 18. Forsman, K.; Glad, T. Constructive algebraic geometry in nonlinear control. Decision and Control, 1990, Proceedings of the 29th IEEE Conference; 1990. p. 2825-2827. 19. Kirk JA, MacGowan GA, Evans C, Smith SH, Warren CM, Mamidi R, Chandra M, Stewart AF, Solaro RJ, Shroff SG. Left ventricular and myocardial function in mice expressing constitutively pseudophosphorylated cardiac troponin I. Circ. Res. 2009; 105:1232–1239. [PubMed: 19850940] 20. Kobayashi T, Solaro RJ. Calcium, thin filaments, and the integrative biology of cardiac contractility. Annu. Rev. Physiol. 2005; 67:39–67. [PubMed: 15709952] 21. Ljung L, Glad T. On global identifiability for arbitrary model parametrizations. Automatica. 1994; 30:265–276. 22. Macgowan GA, Kirk JA, Evans C, Shroff SG. Pressure-calcium relationships in perfused mouse hearts. Am. J. Physiol. Heart Circ. Physiol. 2006; 290:H2614–H2624. [PubMed: 16415077] 23. Peterson JN, Hunter WC, Berman MR. Estimated time course of Ca2+ bound to troponin C during relaxation in isolated cardiac muscle. Am. J. Physiol. 1991; 260:H1013–H1024. [PubMed: 2000960] 24. Pia Saccomani M, Audoly S, D’Angiò L. Parameter identifiability of nonlinear systems: the role of initial conditions. Automatica. 2003; 39:619–632. 25. Rhodes SS, Ropella KM, Audi SH, Camara AK, Kevin LG, Pagel PS, Stowe DF. Cross-bridge kinetics modeled from myoplasmic [Ca2+] and LV pressure at 17 °C and after 37 °C and 17 °C ischemia. Am. J. Physiol. Heart Circ. Physiol. 2003; 284:H1217–H1229. [PubMed: 12531735] 26. Shroff SG, Berger DS, Korcarz C, Lang RM, Marcus RH, Miller DE. Physiological relevance of T-tube model parameters with emphasis on arterial compliances. Am. J. Physiol. 1995; 269:H365– H374. [PubMed: 7631869] 27. Stergiopulos N, Westerhof BE, Westerhof N. Total arterial inertance as the fourth element of the windkessel model. Am. J. Physiol. 1999; 276:H81–H88. [PubMed: 9887020] 28. Toy SM, Melbin J, Noordergraaf A. Reduced models of arterial systems. IEEE Trans. Biomed. Eng. 1985; 32:174–176. [PubMed: 3997173] 29. Walter E, Lecourtier Y. Global approaches to identifiability testing for linear and nonlinear state space models. Math. Comput. Simul. 1982; 24:472–482. 30. Westerhof N, Elzinga G, Sipkema P. An artificial arterial system for pumping hearts. J. Appl. Physiol. 1971; 31:776–781. [PubMed: 5117196]


Kirk et al.

Page 19

Author Manuscript FIGURE 1.

Schematic representation of the example model, showing input, u(t), output (dashed line), y(t), state variables, x1 and x2, and model parameters, k1 and k2.

Author Manuscript Author Manuscript Author Manuscript Cardiovasc Eng Technol. Author manuscript; available in PMC 2015 December 30.

Kirk et al.

Page 20


FIGURE 2.

Author Manuscript

Schematics of systemic arterial circulation models. (a) Electrical analog of the three-element Windkessel model, with two free parameters (C and r or Rp), model input (aortic flow) Qao, and model output (aortic pressure) Pao. (b) Electric analog of the five-element Windkessel model, with four free parameters (r, C1, L, and C2), model input Qao, and model output Pao. (c) Schematic of the T-Tube model with two experimental conditions: (1) seven free parameters (τb, τh, Clb, Clh, Rob, Roh, and Rpb or Rph), model input Pao, and model output Qao and (2) six free parameters (τb, τh, Clb, Clh, Rob, and Roh), model input Pao, and model outputs Qao and Qb. Subscripts b and h denote head and body circulation, respectively and subscript l denotes terminal load.

Author Manuscript Cardiovasc Eng Technol. Author manuscript; available in PMC 2015 December 30.

Kirk et al.

Page 21

Author Manuscript Author Manuscript FIGURE 3.

Author Manuscript

Schematics of the crossbridge cycling-based cardiac muscle contraction models. (a) Threestate model (Xco, Xcn, Xcf) with transition rate constants (kon, koff, f, d, h, g, S), feedback (grey line), model input ([Ca](t)), and model output (Force, F(t)). (b) Four-state model, with transition rate constants (k1, k2, k3, k4, f, g, d, S), feedback (grey lines), model input ([Ca] (t)), and model output (Force, F(t)).

Author Manuscript Cardiovasc Eng Technol. Author manuscript; available in PMC 2015 December 30.

Kirk et al.

Page 22


FIGURE 4.

Pressure waveforms (arbitrary units) for model-generated synthetic data (solid line) and model-based fits (parameter optimization) to dynamic data only (dotted line) or dynamic and steady-state data simultaneously (dashed line).

Author Manuscript Author Manuscript Cardiovasc Eng Technol. Author manuscript; available in PMC 2015 December 30.

Kirk et al.

Page 23

TABLE 1

Author Manuscript

The models that are tested for a priori global identifiability using DAISY, along with the number of state variables, measured outputs, free parameters to estimate, and global identifiability results. # State Var.

# Outputs

# Free parameters

Identifiable?

3-Element

1

1

2

Yes

5-Element

3

1

4

Yes

With one output

2

1

7

No

With two outputs

2

2

6

Yes

Stiffness model

2

1

4

Yes

Without feedback

3

1

8

Yes

With feedback

3

1

9

Yes

Without feedback

3

1

8

Yes

With feedback

3

1

10

Yes

Model Windkessel model

T-tube model

3-State model

Author Manuscript

4-State model


Kirk et al.

Page 24

TABLE 2

Author Manuscript

Numerical parameter estimation analysis for the four-state model of cardiac muscle contraction.

Author Manuscript

Parameter

Actual parameter value used to create synthetic force data

k2

25

3.42 ± 144

24.6 ± 0.044

k3

84

93.3 ± 13.3

86.2 ± 0.176

k4

84

80.0 ± 55.9

82.9 ± 0.147

g

205.3

201 ± 76.3

206 ± 0.349

d

44.6

42.5 ± 4.07

44.0 ± 0.067

α1

2

1.50 ± 15.9

2.04 ± 0.003

β1

3

2.63 ± 15.0

3.05 ± 0.003

αf

0.15

0.130 ± 1.18

0.150 ± 0.001

βf

0.7

0.839 ± 4.78

0.705 ± 0.001

N/A

>0.99

>0.99

R

2

Estimated parameter fit to synthetic force data (dynamic)

Estimated parameter fit to synthetic force data (dynamic + steady-state)

Parameters values used to generate a force transient corresponding to a typical calcium transient (dynamic data). Additionally, these parameter values were used to generate a set of steady-state data, in which a constant calcium concentration was applied to the model, and force was recorded after the steady-state value was achieved. This was done for a range of constant calcium levels (from zero to saturating conditions). Numerical parameter estimation was performed on either the dynamic data alone (third column, estimated parameter value ±95% confidence interval) or both dynamic and steady-state data simultaneously (fourth column, estimated parameter value ±95% confidence interval). The gray boxes indicate those parameters that were estimated relatively accurately (i.e., with small confidence interval) when fit to only the dynamic synthetic data.


Kirk et al.

Page 25

TABLE 3

Author Manuscript

Sensitivities (mean ± SEM) of selected calculated indices to changes in parameter values. Calculated index (steady-state data) F max

pCa50

k2

−0.04 ± 0.04

−0.03 ± 0.01

−0.46 ± 0.19

−0.05 ± 0

−0.12 ± 0.15

0.03 ± 0.04

k3

−0.05 ± 0.03

−0.11 ± 0.03

−0.42 ± 0.17

−0.94 ± 0.38

−0.42 ± 0.15

−0.09 ± 0.11

Parameter

nH

Calculated index (dynamic data) F dev

T rise

T relax

k4

0.03 ± 0

0.04 ± 0

0.46 ± 0.09

0.29 ± 0.02

0.08 ± 0.29

0.01 ± 0.02

g

−1.02 ± 0.29

−0.11 ± 0.02

−1.65 ± 0.84

−0.83 ± 0.3

−0.12 ± 0.15

0.03 ± 0.04

α1

0.03 ± 0.01

0.08 ± 0.01

1.01 ± 0.19

0.05 ± 0

0.12 ± 0.15

−0.03 ± 0.04

Author Manuscript

β1

0±0

0.03 ± 0

−0.77 ± 0.45

1.01 ± 0

0.12 ± 0.15

0±0

αf

0.58 ± 0.01

0.04 ± 0.01

1.82 ± 1.11

0.04 ± 0

0.12 ± 0.15

0±0

βf

0.41 ± 0.08

0.09 ± 0.02

0.21 ± 0.08

1.01 ± 0.01

0.12 ± 0.15

−0.03 ± 0.04

d

−0.01 ± 0

−0.04 ± 0.01

−0.23 ± 0.07

−0.3 ± 0.08

−0.23 ± 0.15

−0.44 ± 0.25

There are three indices calculated from the steady-state force output (Fmax: maximal calcium-activated force, pCa50: calcium which generates 50% Fmax, and nH: Hill coefficient, or steepness of the sigmoidal curve) and three indices calculated from the dynamic force output (Fdev: peak developed force, Trise: rise time, Trelax: relaxation time). Sensitivity values were calculated by altering parameters over the range of 50 to 150% of the baseline values (see text).


Structural Identifiability of Dynamic Systems Biology Models.

Parameterisation affects identifiability of population models.

On the Identifiability of Transmission Dynamic Models for Infectious Diseases.

Robust global identifiability theory using potentials--Application to compartmental models.

Local identifiability and sensitivity analysis of neuromuscular blockade and depth of hypnosis models.

On linear models and parameter identifiability in experimental biological systems.

Practical identifiability of growth and substrate consumption models.

A review of a priori regression models for warfarin maintenance dose prediction.

Determination of parameter identifiability in nonlinear biophysical models: A Bayesian approach.

Mediterranean Diet and Cardiovascular Disease: A Critical Evaluation of A Priori Dietary Indexes.

Structural and practical identifiability analysis of S-system.

Comparison of approaches for parameter identifiability analysis of biological systems.

Identifiability of models for time-resolved fluorescence with underlying distributions of rate constants.

Statistical identifiability and convergence evaluation for nonlinear pharmacokinetic models with particle swarm optimization.

Structural identifiability of cyclic graphical models of biological networks with latent variables.

Clustering reveals limits of parameter identifiability in multi-parameter models of biochemical dynamics.

Identifiability and parameter estimation.

Ross's a priori pathometry - a perspective.

Predictive a priori pressure-dependent kinetics.

A computer program for the analysis of controllability, observability and structural identifiability of biological compartmental systems.

A priori estimation of memory effects in reduced-order models of nonlinear systems using the Mori-Zwanzig formalism.

Template based rotation: a method for functional connectivity analysis with a priori templates.

Parameter estimation: local identifiability of parameters.

Optimizing cardiovascular benefits of exercise: a review of rodent models.