## Abstract

The study assesses the strength of the causal relation along baroreflex (BR) in humans during an incremental postural challenge soliciting the BR. Both cardiac BR (cBR) and sympathetic BR (sBR) were characterized via BR sequence approaches from spontaneous fluctuations of heart period (HP), systolic arterial pressure (SAP), diastolic arterial pressure (DAP) and muscle sympathetic nerve activity (MSNA). A model-based transfer entropy method was applied to quantify the strength of the coupling from SAP to HP and from DAP to MSNA. The confounding influences of respiration were accounted for. Twelve young healthy subjects (20–36 years, nine females) were sequentially tilted at 0°, 20°, 30° and 40°. We found that (i) the strength of the causal relation along the cBR increases with tilt table inclination, while that along the sBR is unrelated to it; (ii) the strength of the causal coupling is unrelated to the gain of the relation; (iii) transfer entropy indexes are significantly and positively associated with simplified causality indexes derived from BR sequence analysis. The study proves that causality indexes are complementary to traditional characterization of the BR and suggests that simple markers derived from BR sequence analysis might be fruitfully exploited to estimate causality along the BR.

This article is part of the themed issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

## 1. Introduction

Baroreflex (BR) plays an important role in the homeostatic regulation of cardiovascular variables in humans [1]. Indeed, in hypovolaemic situations that might lead to an arterial pressure (AP) drop, such as standing, the activation of this reflex triggers tachycardic responses of the heart and vasoconstriction at the periphery to maintain appropriate AP values.

The cardiac BR (cBR) evokes heart period (HP) responses to AP variations [2]. The cBR can be characterized from spontaneous variability of HP and systolic AP (SAP) by extracting parallel sequences of consecutive HP and SAP variations (i.e. positive variations of SAP accompanied by positive changes of HP or *vice versa*), usually referred to as cBR sequences [3]. The cBR sensitivity (cBRS), estimated as the average slope of the linear regression computed in the (SAP,HP) plane over the cBR sequences, the percentage of cBR sequences with respect to the overall amount of sequences and the fraction of cBR sequences with respect to the overall amount of positive and negative SAP ramps are traditional descriptors of cBR [3–5]. The strength of the causal relation from SAP to HP estimated via a Granger–Wiener causality approach is a more recently proposed marker of cBR regulation [6,7].

The sympathetic BR (sBR) governs the responses of sympathetic activity to AP changes [8]. The sBR can be typified in humans from spontaneous variability of diastolic AP (DAP) and sympathetic discharge. In contrast to cBR, a minimally invasive set-up for acquiring muscle sympathetic nerve activity (MSNA), usually from the peroneal nerve, is necessary [8]. A suitable definition of sBR sequences as antiparallel successions of consecutive MSNA burst rate and SAP variations (i.e. positive variations of DAP accompanied by negative changes of MSNA bust rate or *vice versa*) allows the calculation of the sBR sensitivity (sBRS) as the average slope of the linear regression computed in the (DAP,MSNA) plane over the sBR sequences and the percentage of sBR sequences with respect to the overall amount of sequences [9]. Similarly to the cBR characterization, the level of sBR solicitation might be estimated through the computation of the percentage of sBR sequences with respect to the overall amount of positive and negative DAP ramps and the strength of the causal relation from DAP to MSNA burst rate variability assessed via a Granger–Wiener causality approach.

Respiration is a confounding factor for both HP-SAP and MSNA-DAP dynamical relations. Indeed, HP variations at the respiratory rate, usually referred to as respiratory sinus arrhythmia [10], might occur even in absence of synchronous modifications of SAP as a result of the modulation of the motoneuron vagal responsiveness by respiratory centres [11]. As a matter of fact, disambiguating direct influences of respiration on HP is mandatory for a reliable assessment of cBR in both animals [12] and humans [13]. Similarly, MSNA burst rate modulations at the respiratory rate might be present in absence of respiratory-related DAP fluctuations. Indeed, the magnitude of the respiratory-related MSNA burst rate modulations is usually considerable [14], while synchronous modifications of DAP are negligible [15]. In addition, it was observed that a MSNA burst is more likely to be found during expiration than during inspiration even in absence of respiratory-related modifications of DAP [16] and this phenomenon is potentiated during deep, low frequency breathing regardless of the DAP levels [17].

The aim of this study is to compute the strength of both cBR and sBR via a Granger–Wiener causality approach during an incremental postural challenge and to correlate it with more standard descriptors of cBR and sBR physiology. A model-based causality approach assessing transfer entropy was exploited to quantify the strength of the coupling from SAP to HP and from DAP to MSNA while disambiguating the direct effect of respiration on both relations [7] and BR sequence methods were used to typify cBR and sBR [3,4,9]. An incremental postural challenge of limited intensity, performed via graded head-up tilt with inclination up to 40° [18], was conducted because this stressor is well-known to decrease respiratory sinus arrhythmia [19,20], increase tonic sympathetic activity and its modulation [21–23] and decrease cBRS and sBRS in proportion to the stimulus [9,21].

## 2. Methods

### (a) Model-based assessment of the strength of causal relation from a source to a target in presence of a confounding factor

We exploited the model-based approach described in [7] for assessing the strength of the causal link from a source signal *x _{i}* = {

*x*(

_{i}*k*),

*k*= 1, … ,

*N*} to a target one

*y*= {

*y*(

*k*),

*k*= 1, … ,

*N*} in presence of a confounding factor

*x*= {

_{c}*x*(

_{c}*k*),

*k*= 1, … ,

*N*}, where

*k*is the progressive counter and

*N*is the series length. Given the full universe of knowledge

*Ω*= {

*y*,

*x*,

_{i}*x*}, the strength of the relation from

_{c}*x*to

_{i}*y*in

*Ω*is assessed via the conditional joint transfer entropy (CJTE) from

*x*and

_{i}*x*to

_{c}*y*given

*x*evaluating the portion of information transferred from

_{c}*x*and

_{i}*x*to

_{c}*y*that can be genuinely attributed to

*x*[24–30]. Under the hypothesis of Gaussianity, this quantity can be computed by comparing the variances of the prediction error in the restricted

_{i}*Ω*obtained from the full

*Ω*by excluding the presumed cause

*x*(i.e.

_{i}*Ω*∖{

*x*} = {

_{i}*y*,

*x*}) and in the full

_{c}*Ω*, and respectively, via the log-ratio [7,26] as 2.1

The variances of the prediction error were computed after fitting *y* in *Ω*∖{*x _{i}*} and

*Ω*via a multivariate linear regression model [12,31]. All model coefficients were identified via traditional least-squares approach and Cholesky decomposition method [32]. All regressions exhibited the same number of coefficients referred to as model order. The model order was optimized in the range from 8 to 14 according to the Akaike figure of merit for multivariate processes [33] over the most complex model structure (i.e. the model of

*y*in

*Ω*). The model order was maintained in

*Ω*∖{

*x*} but the coefficients of the models were estimated again. After estimating the model parameters, the best one-step-head prediction of the current value of

_{i}*y*can be computed and the prediction error can be calculated as the difference between the original value and its best prediction [13]. The calculation of the variances of the prediction error in

*Ω*∖{

*x*} and

_{i}*Ω*allows the computation of (2.1).

### (b) Traditional characterization of cBR via the sequence method

cBRS was estimated according to cBR sequence method [3] as implemented in [12]. This technique relies on the search for spontaneous HP and SAP sequences characterized by the contemporaneous increase or decrease of HP and SAP. Sequences were composed of four consecutive HPs and SAPs (i.e. three variations) that met the following prerequisites [5]: (i) the absolute value of the total HP variation was larger than 5 ms; (ii) the absolute value of the total SAP variation was larger than 1 mm Hg; (iii) the linear correlation coefficient, r_{HP‐SAP}, computed in the [SAP(*i*), HP(*i* + *τ*_{HP‐SAP})] plane over the correspondent HP and SAP sequences, was larger than 0.85, where *τ*_{HP‐SAP} represents the lag between HP and SAP. *τ*_{HP‐SAP} was set to 0 to pick up the fast vagal arm of the cBR capable of modifying HP in reaction of SAP changes occurring within the same HP [34,35]. The slope of the regression line was calculated for each cBR sequence and subsequently averaged over all sequences. This average was taken as an estimate of cBRS and expressed in ms · mm Hg^{−1}. The percentage of cBR sequences (cBRS%) was computed as the number of cBR sequences multiplied by 100 and divided by the total number of sequences. The cBR effectiveness index (cBEI) was defined as the number of cBR sequences divided by the total number of positive and negative SAP ramps [4]. Both cBRS% and cBEI are dimensionless.

### (c) Characterization of sBR via the sequence method

sBRS was estimated according to the sBR sequence method [9]. The sBR sequence technique relies on the search for spontaneous MSNA burst rate and DAP sequences characterized by the contemporaneous increase of MSNA burst rate and decrease of DAP or *vice versa*. The sequences were composed of four consecutive MSNA burst rate and DAP values (i.e. three variations) that met the following prerequisites [9]: (i) the absolute value of the total MSNA burst rate variation was larger than 0 bursts · s^{−1}; (ii) the absolute value of the total DAP variation was larger than 1 mm Hg; (iii) the linear correlation coefficient, *r*_{MSNA‐DAP}, computed in the [DAP(*i*), MSNA(*i* + *τ*_{MSNA‐DAP})] plane over the correspondent MSNA burst rate and DAP sequences, was larger than 0.85, where *τ*_{MSNA‐DAP} represents the lag between MSNA burst rate and DAP. *τ*_{MSNA‐DAP} was set to 0 [9] fully in agreement with the typical latency of sBR [36,37] given the definition of MSNA burst rate variability [23]. The slope of the regression line was calculated for each sBR sequence and subsequently averaged over all sequences. This average was taken as an estimate of sBRS and expressed in bursts s^{−1 }mm Hg^{−1}. The percentage of sBR sequences (sBRS%) and the sBR effectiveness index (sBEI) were calculated respectively as the number of sBR sequences multiplied by 100 and divided by the total number of sequences and as the number of sBR sequences divided by the total number of positive and negative DAP ramps. Both sBRS% and sBEI are dimensionless.

## 3. Experimental protocol, series extraction and analysis of data

### (a) Experimental protocol

The study comprised 12 healthy subjects (age from 20 to 36 years, median = 22.5 years; body mass index from 18.6 to 28.4 kg m^{−2}, median = 24.2 kg m^{−2}; nine females). The study protocol was approved by the Alfred Hospital Ethics Review Committee (n. 144/06) and conformed to the relevant guidelines of the National Health and Medical Research Council of Australia and to the principles of the Declaration of Helsinki for medical research involving humans. All subjects provided written informed consent. The criteria for subjects' selection and the description of the head-up tilt protocol were provided in [18]. Briefly, the surface electrocardiogram (ECG), invasive AP from the radial artery, respiratory movements (RMs) via a piezoelectric device and multiunit MSNA in postganglionic fibres distributed to the skeletal muscle vasculature from the common peroneal nerve were simultaneously recorded during head-up tilt with tilt table inclination of 0°, 20°, 30° and 40° (T0, T20, T30 and T40 respectively) for 10 min at each angle. The head-up tilt test started from the horizontal position and was incremental with respect to the previous tilt table inclination. The raw MSNA signal was band-pass filtered (700–2000 Hz), amplified, rectified and integrated (time constant of 0.1 s) to obtain integrated MSNA. ECG, AP, RM and integrated MSNA were digitized at 1000 Hz using a PowerLab system (model ML785/8SP, ADInstruments, Castle Hill, NSW, Australia) and stored for off-line analysis. Out of all recordings obtained from the 12 subjects, signals were of insufficient quality or the tilt protocol was not completed, respectively, for one subject during T30 and T40. At difference with [9] recordings during T60 were not used in this study given the inability of finding stationary RM traces in this specific experimental condition. The same situation was found in four, two, one and one subjects during T0, T20, T30 and T40. Full access to this database is available free of charge by contacting the corresponding author.

### (b) Extraction of the beat-to-beat variability series

MSNA bursts were detected from the integrated MSNA signal by using an adaptive thresholding method and searching for the occurrence of the MSNA burst from 0.9 to 1.7 s from the current R-wave peak according to the latency of the sBR [36,37]. Then, the MSNA bursts associated to an R-wave peak were counted over a moving time window of 5 s that was advanced in steps of 1 ms along the MSNA recording. The resulting step-wise burst-count MSNA signal was divided by the window length and low-pass filtered with a finite impulse response filter with a cut-off frequency of 0.5 Hz. This signal represents the modulation of the MSNA burst rate, its values are expressed in bursts s^{−1} and has the same time resolution as the original integrated MSNA signal [23]. The HP, SAP, DAP, RM and MSNA measures were taken on a beat-to-beat basis as follows: (i) the current HP was measured as the temporal distance between two consecutive R-wave peaks detected on the ECG; (ii) the *k*th SAP value was taken as the maximum of AP in the *k*th HP; (iii) the *k*th DAP was computed as the minimum of the AP preceding the *k*th SAP; (iv) the *k*th RM measure was taken by sampling the RM signal at the first R-wave delimiting the *k*th HP; (v) the *k*th MSNA burst rate value was obtained by sampling the signal representing the modulation of the MSNA burst rate at the first R-wave delimiting the *k*th HP [23]. HP, SAP, DAP, RM and MSNA burst rate beat-to-beat variability series of *N* = 256 consecutive values were randomly selected within each experimental condition and linearly detrended. If non-stationarities, such as very slow drifting of the mean or sudden changes of the variance, were detected using a restricted weak stationarity test [38], the random selection was repeated.

### (c) Computation of causality markers

In this work the directed information transfer indexes along the cBR were computed in *Ω* = {HP, SAP, RM} and those along the sBR were calculated in *Ω* = {MSNA, DAP, RM}. The AP series (i.e. SAP and DAP) was considered the source signal given that AP is the sensed variable of the BR. RM is the conditioning variable because it perturbs both HP-SAP and MSNA-DAP dynamical relations. HP and MSNA are the target variables of cBR and sBR respectively. In the computation of the delay from SAP and RM to HP was set to 0 to account for a fast vagal action of cBR and cardiopulmonary reflexes [13,35]. In the computation of the delay from DAP and RM to MSNA burst rate was set to 0 [9] because the duration of the running window used for the calculation of the MSNA burst rate (i.e. 5 s) makes this value fully compatible with the latency of sBR and disturbing action of respiratory activity [36,37]. We computed and as well as cBRS, cBRS%, cBEI and sBRS, sBRS%, sBEI.

### (d) Statistical analysis

The linear association of the descriptors of the cBR and sBR with the magnitude of the stimulus was tested by checking whether the slope of the linear regression of the parameters on the component of the gravity parallel to the tilt table was significantly different from 0. The same analysis was exploited to check the association between the strength of the causal relation and the correspondent traditional descriptors of the reflex after pooling together all subjects regardless of the experimental condition. Pearson product moment correlation *r* and type I error probability *p* were reported. Statistical analysis was carried out using a commercial statistical program (Sigmaplot, Systat Software, Inc, Chicago, IL, USA, v.11.0). A *p* < 0.05 was considered statistically significant.

## 4. Results

Traditional markers of cBR (i.e. cBRS and cBRS%) were significantly linearly related to tilt table angles: the correlation coefficient of cBRS was negative (i.e. *r* = −0.474, *p* = 2.67 × 10^{−3}), while that of cBRS% was positive (i.e. *r* = 0.378, *p* = 1.94 × 10^{−2}). The descriptors of sBR (i.e. sBRS and sBRS%) were significantly linearly related to magnitude of the postural challenge as well, but both the indexes were positively related to tilt table angles (i.e. *r* = 0.343, *p* = 3.49 × 10^{−2} and *r* = 0.457, *p* = 3.90 × 10^{−3} respectively).

Figure 1 shows the individual values (open circles) of (figure 1*a*) and (figure 1*b*) as a function of the sine of the tilt table inclination. The regression line of the variable on the sine of the tilt table angles (solid line) and its 95% confidence interval (dotted lines) are plotted as well, if the slope of the regression line is significantly larger than 0. A significant positive association of with the magnitude of the postural challenge was found (figure 1*a*, *r* = 0.355, *p* = 2.89 × 10^{−2}), while was unrelated to tilt table inclination (figure 1*b*, *r* = 0.170, *p* = 3.09 × 10^{−1}).

Figure 2 explores the presence of linear associations between directed information flow indexes and the gain of the relation. The panels show the pairs (open circles) in the (cBRS, ) plane (figure 2*a*) and in the (sBRS, ) plane (figure 2*b*), estimated from each subject regardless of the experimental condition. The regression lines (solid line) and their 95% confidence interval (dotted lines) are plotted as well if the slope of the regression line is significantly larger than 0. No significant association between causality indexes and gains of the reflex was detected (*r* = −0.239, *p* = 1.48 × 10^{−1} and *r* = −0.0954, *p* = 5.69 × 10^{−1} respectively).

Figure 3 explores the presence of linear associations between indexes of the strength of the causal link and the percentage of BR sequences (figure 3*a,b*) and between the strength of the causal link and the fraction of ramps of the input (i.e. SAP or DAP) inducing a variation of the output compatible with a working BR arm (figure 3*c*,*d*). Accordingly, the panels show the pairs (open circles) in the planes (cBRS%, ) (figure 3*a*), (sBRS%, ) (figure 3*b*), (cBEI,) (figure 3*c*) and (sBEI,) (figure 3*d*) as estimated from each subject regardless of the experimental condition. The regression lines (solid line) and their 95% CI (dotted lines) are plotted as well if the slope of the regression line is significantly larger than 0. Both and were significantly and positively correlated with their correspondent cBEI (figure 3*c*, *r* = 0.464, *p* = 3.35 × 10^{−3}) and sBEI (figure 3*d*, *r* = 0.611, *p* = 4.59 × 10^{−5}). Conversely, while was significantly associated with sBRS% (figure 3*b*, *r* = 0.520, *p* = 8.28 × 10^{−4}), was unrelated to cBRS% (figure 3*a*, *r* = 0.262, *p* = 1.12 × 10^{−1}).

## 5. Discussion

The main findings of this study can be summarized as follows: (i) even in response to a postural stimulus of limited intensity the strength of the causal relation from SAP to HP along the cBR increases as a function of the tilt table inclination; (ii) the strength of the causal link from DAP to MSNA is unrelated to the magnitude of the postural challenge; (iii) the strength of the causal coupling is unrelated to the gain of the dynamical relation and this conclusion holds for both cBR and sBR; (iv) causality indexes along cBR and sBR are significantly and positively associated with cBEI and sBEI respectively.

### (a) The information transfer from SAP to HP gradually increases with the magnitude of the postural challenge even in presence of a stimulus of a limited intensity

It is well-known that the strength of the causal link from SAP to HP along the cBR augments with the magnitude of the postural stimulus [6,7,13,39,40], while cBRS gradually decreases [9,21,42]. The former is an indication of the increased solicitation of cBR with the magnitude of the postural stressor, observable even at the level of cBRS% [9], while the latter is the result of the vagal withdrawal and sympathetic activation decreasing respiratory sinus arrhythmia and augmenting the amplitude of the SAP fluctuations in proportion to the stimulus [19–23]. Indeed, the same conclusion was reached by both a model-free bivariate approach disregarding RM [6,39] as well as via a model-based trivariate approach accounting for the RM as a latent confounder affecting both SAP and HP [7,13,40]. In this study the information genuinely transferred from SAP to HP along the cBR, taken as a marker of the intensity of the directed action from SAP to HP [24–30], increased with tilt table angle even in presence of a stimulus of limited intensity. This finding suggests that even limited orthostatic challenges can solicit cBR and the threshold for its activation, if present, is very small.

### (b) The information transfer from DAP to MSNA burst rate remains stable with the magnitude of the postural challenge

Results of the present study confirm that, even in presence of a stimulus of limited intensity, sBRS moves towards 0 and sBRS% increases with tilt table angle [9]. Despite the progressive modification of sBRS and sBRS% the strength of the causal relation from DAP to MSNA burst rate along the sBR remains steady as a function of the tilt table inclination. Given that the dependency of the strength of the causal link from DAP to MSNA burst rate is expected as a sign of the greater solicitation of sBR with the stimulus, we suggest that, in this specific experimental protocol exploiting a stimulus of limited intensity, the relevance of the stressor is not sufficient to reveal the progression of the genuine information transfer from DAP to MSNA burst rate. We conjecture that, given the low intensity of the postural challenge, MSNA burst rate might be predominantly under central commands [41] operating independently of DAP levels and to a lesser extent reflexively regulated via sBR, thus masking the expected positive relation of the strength of the directed information flow with tilt table inclination. In addition, taking together the gradual increase of the strength of the causal relation from SAP to HP and the steady behaviour of that from DAP to MSNA burst rate, we conclude that cBR is more responsive in regulating AP following postural challenge than sBR. Different thresholds for activation of cBR and sBR, with cBR becoming operative with more limited stressors than sBR, might explain this behaviour as well.

### (c) The strength of the causal coupling is unrelated with the gain of the dynamical relation

From a methodological standpoint, the magnitude of the input–output dynamical relation provides information completely different from the strength of the causal relation from the input to the output [43–45]. The magnitude of the input–output dynamical relation, usually computed under the hypothesis of linearity and time invariance as the modulus of the transfer function from the input to the output, provides the ratio of the amplitude of the output oscillation to that of the input oscillation at the same frequency as a function of the frequency [46]. As a consequence of its definition, the dimension of the magnitude of the dynamical relation is the ratio of the dimensions of the input to that of the output. The strength of the causal relation from the input to the output is naturally computed via metrics assessing the association between the input and the output in a given temporal direction, such as predictability improvement and transfer entropy, even though techniques in the frequency domain are frequently used [30]. As a consequence of its definition, the dimension of the strength of the causal link from the input to the output is dimensionless and usually bounded. Given this premise, it is not surprising to find measures of the strength of the causal relation along cBR and sBR unrelated to cBRS and sBRS respectively. This finding suggests that causality analysis provides non-redundant information compared to more traditional analyses involving the assessment of the gain of the input–output relation and should be used in association with them to exploit its complementary view on the data.

### (d) cBR and sBR sequences analyses provide markers of causality

Nollo *et al.* [47] suggested that cBEI [4] is an index of the strength of the causal link along the cBR. Indeed, cBEI, by selecting among the set of positive and negative SAP ramps (i.e. the effective input for the cBR) exclusively those producing an effective parallel modification of HP, is capable of quantifying the HP-SAP association in the temporal direction of the cBR. In agreement with this observation, the genuine information transfer from SAP to HP is significantly and positively correlated to cBEI. Remarkably, this result holds even for sBR, being the genuine information transfer from DAP to MSNA burst rate significantly and positively associated with sBEI. This result suggests that simplified indexes of causality might be computed without needing the identification of a model [26,30,48], or demanding calculus, as in model-free approaches [24,25,27,29,30]. The suitability of this simplified approach might be limited to situations in which the input–output relation can produce parallel or antiparallel responses of the effect signal to a given variation of the presumed cause, like in the case of cBR and sBR respectively. Further tests are needed to understand whether this basic approach holds in pathological populations characterized by an impaired or missing cBR and sBR, especially in presence of unaltered feed-forward mechanisms assuring association between the two variables in the reverse temporal direction and in presence of more complex responses of the target variables to increases or decreases of the presumed cause. In this study model-based transfer entropy indexes are computed by conditioning out respiration, while simplified markers of causality based on BR sequence analysis ignore its influence. Given that respiration is a latent confounder for both cBR and sBR directed links [7,13,40], the high correlation between the two types of causality indexes might be surprising. A possible explanation of this finding is that responses of HP and MSNA burst rate to positive and negative AP ramps occur regardless of the respiratory phases, thus making irrelevant to account for respiration in this specific protocol.

## 6. Conclusion

The study assesses the strength of the causal relation along cBR and sBR via a model-based transfer entropy approach during an incremental orthostatic challenge of limited intensity. It suggests a potential difference between the two arms of the BR, given that cBR is activated in proportion with the magnitude of the stressor even in presence of a stimulus of limited magnitude, while the hypothesized progressive activation of sBR is not visible even though sBRS moves gradually towards 0. Since indexes of causality are unrelated to the gain of the reflex, the study proves experimentally their complementary value in the field of cardiovascular control analysis based on spontaneous variability, even though further studies are necessary to clarify the practical usefulness of causality markers. Moreover, since simplified indexes of causality based on BR sequence analysis are highly correlated with more sophisticated model-based information flow markers, this study suggests the possibility to perform causality analysis as a by-product of the more standard characterization of the BR via sequence analysis. We advocate applications to pathological population to validate further this conjecture and to better define the range of the practical applicability of these simplified causality indexes.

## Authors' contributions

A.P. conceived the study; A.P. drafted the manuscript; A.P. and A.M. analysed the data; M.E. and E.L. performed the experiments; A.P., A.M., V.B., B.D.M., M.E., E.L. and M.B. interpreted the data; A.P., A.M., V.B., B.D.M., M.E., E.L. and M.B. critically revised the manuscript; A.P., A.M., V.B., B.D.M., M.E., E.L. and M.B. approved the final version of the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

No financial support was received for this research.

## Footnotes

One contribution of 14 to a theme issue ‘Mathematical methods in medicine: neuroscience, cardiology and pathology’.

- Accepted December 9, 2016.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.