- Open Access
Age and gender dependent heart rate circadian model development and performance verification on the proarrhythmic drug case study
Theoretical Biology and Medical Modellingvolume 10, Article number: 7 (2013)
There are two main reasons for drug withdrawals at the various levels of the development path – hepatic and cardiac toxicity. The latter one is mainly connected with the proarrhythmic potency and according to the present practice is supposed to be recognized at the pre-clinical (in vitro and animal in vivo) or clinical level (human in vivo studies). There are, although, some limitations to all the above mentioned methods which have led to novel in vitro – in vivo extrapolation methods being introduced. With the use of in silico implemented mathematical and statistical modelling it is possible to translate the in vitro findings into the human in vivo situation at the population level. Human physiology is influenced by many parameters and one of them which needs to be properly accounted for is a heart rate which follows the circadian rhythm. We described such phenomenon statistically which enabled the improved assessment of the drug proarrhythmic potency.
A publicly available data set describing the circadian changes of the heart rate of 18 healthy subjects, 5 males (average age 36, range 26–45) and 13 females (average age 34, range 20–50) was used for the heart rate model development. External validation was done with the use of a clinical research database containing heart rate measurements derived from 67 healthy subjects, 34 males and 33 females (average age 33, range 17–72). The developed heart rate model was then incorporated into the ToxComp platform to simulate the impact of circadian variation in the heart rate on QTc interval. The usability of the combined models was assessed with moxifloxacin (MOXI) as a model drug.
The developed heart rate model fitted well, both to the training data set (RMSE = 128 ms and MAPE = 12.3%) and the validation data set (RMSE = 165 ms and MAPE = 17.1%). Simulations performed at the population level proved that the combination of the IVIVE platform and the population variability description allows for the precise prediction of the circadian variation of drugs proarrhythmic effect.
It can be concluded that a flexible and practically useful model describing the heart rate circadian variation has been developed and its performance was verified.
The withdrawal of several marketed drugs in the last decade, extensive re-labeling and a high attrition rate, especially in the late stage of the drug development process, due to the QT prolonging liability or TdP occurrence, have all established that proarrhythmia assessment is a primary focus of regulatory bodies and the pharmaceutical industry. Thus careful evaluation of the proarrhythmic potential of an investigated compound is an integral element of the safety profile required for the approval of new drugs. There are several available models utilized for the assessment of drugs proarrhythmic potency at various levels of drug development. The most commonly utilized and probably most reasonable pathways at the non-clinical stage of development include in vitro models where the human ionic channels are heterologously expressed in non-human (i.e. XO, CHO) and human (i.e. HEK) cells, in vivo/ex vivo animal studies where suggested species are used for the drug influence evaluation at the cardiac level . Novel techniques include the human stem cells derived cardiomyocytes which, however, still lack the standard methodology . At the early clinical development stage the currently best available, yet costly and still not perfect thorough QT/QTc clinical studies (TQTs)  are introduced to assess the QT prolongation by the candidate drug as compared with placebo and positive control. Currently there is a translational gap in the quantitative extrapolation of the in vitro or animal studies’ outcomes to the human situation. There are some rules of thumb  and decision trees  available to qualitatively (yes/no) predict QT prolongation potential in humans from in vitro assays, particularly hERG channel inhibition assays. Moreover, there are no methods available to combine in vitro assay information from multiple ion channels (INa, ICa, IKr, IKs) with human physiological and demographic data to predict QT prolongation in humans. Physiologically based in silico methods have a potential to bridge the translation gap between the in vitro electrophysiological data gained from in vitro/ex-vivo studies and human QT liability. We have recently developed a so called “bottom-up” in silico platform (ToxComp) which combines a physiologically based electrophysiological model of human left ventricular cardiomyocytes and a database of human physiological, genotypic and demographic data enabling the prediction of the QT prolongation in humans based on the in vitro data [6, 7]. ToxComp has been previously used to predict QT prolongation liability of various compounds [8, 9]. Until now ToxComp has predicted the absolute QT prolongation liability of a compound assuming that the physiological parameter values are constant during the simulation, which limited its ability to simulate clinical scenarios. With the addition of the circadian variation model in ToxComp, the platform could potentially be employed in simulating clinical scenarios before conducting clinical studies with an aim to optimize the designs and reduce the chances of the failure of these costly and time intensive clinical studies.
A circadian rhythm is a biological process that displays an endogenous, entrainable oscillation of roughly 24 hours . In humans there are multiple physiological processes which follow the circadian clock, including the heart rate (HR) . In our study context, however, the circadian rhythm of the heart rate is one of the most important, since it may both influence the QT length and interfere with the drug effect . The mathematical and statistical modelling of the circadian variation of the heart rate has been previously studied. Nakagawa and colleagues reported results of a study which involved 44 (but 20 were included into the analysis phase) healthy individuals. Subjects showed a significant circadian heart rate rhythm and based on the collected data the authors proposed a cosinor model to describe such phenomenon, although they did not differentiate between male and females which limits the practical use of the derived model . Similar models, carrying identical limitations were proposed independently by Massin and colleagues and Li and colleagues [14, 15]. 57 healthy children and 115 healthy, non-smoking adults were included respectively and the 24-hours ECG measurements were used to derive the regression models. Probably the most detailed analysis of the heart rate variation was presented by Bonnemeier et al. . A large group of 166 healthy individuals (81 females, 85 males) characterized by a wide age range (20–70 years) was studied. The authors analysed differences of heart rate circadian profiles for hourly aggregated measurements in a groups stratified by age decades, separately for female and male subjects. The authors did not, however, model their data. There is also a large number of publications where circadian heart rate variation is discussed among various subpopulations (athletes, truck drivers, welders) and diseases (myotonic dystrophy, angor patients) [17–20]. Although to our best knowledge none of them proposed a model flexible enough and described in enough detail to be directly applicable for the generation of a virtual human population where heart rate is an attribute specific for every individual, and this was the main reason we decided to develop a new model.
Aims of the study
We aimed to develop and validate a model describing the gender and age dependent circadian rhythm of the heart rate in European Caucasian healthy subjects. Additional aims included quality and usability assessment of the model when incorporated into the ToxComp platform. Its usability, in combination with other parameters describing demographic (age, gender), physiological (i.e. cardiomyocytes characteristics, plasma ions concentration) and genetic variability in a population, was tested in the virtual clinical trial.
Materials and methods
PhysioNet data set description
The analyzed data set was obtained from the PhysioBank, which is an archive of digitized sets of data reflecting physiological signals. The data warehouse contains over 50 various freely available databases, and for modelling purposes the MIT-BIH Normal Sinus Rhythm Database was used . There were a total of 18 subjects, 5 males (average age 36, range 26–45) and 13 females (average age 34, range 20–50). For each subject up to 24 hours of RR recordings were available. Subjects, on average, had 94,440 individual RR measurements (range 73,300-115,900). In order to decrease the significant computational burden it was necessary to reduce the number of RR measurements per subject. RR averaging in one-hour, half-hour and 15-minute ‘time windows’ was found inadequate due to its variability reduction property. Stable results were obtained after sampling RR measurements every 1 minute.
Validation data set
Model validation was performed with the use of a completely independent data set. Data was derived from Cracow’s clinical research database (1st Department of Cardiology and Hypertension, Jagiellonian University Medical College). There were a total of 67 healthy subjects in the validation data set, 34 males and 33 females (average age 33, range 17–72). Validated oscillometric SpaceLabs 90207 monitors (Redmond, WA, USA) fitted with the appropriate cuff size were programmed to obtain readings every 15 minutes from 08.00 to 22.00 and every 30 minutes from 22.00 to 08.00. Each reading included systolic/diastolic blood pressure, mean arterial pressure and heart rate. Subjects, on average, had 71 heart rate measurements (range 35–99). RR interval length was calculated based on the heart rate measurement results.
Model usability testing
To test the usability of the implemented circadian rhythm model within the ToxComp platform, we have chosen moxifloxacin (MOXI), the most commonly used positive control in TQT studies as the model compound. The ToxComp system was used to simulate the drugs triggered ECG modification. The ToxComp platform combines a physiologically based electrophysiological model of human left ventricular cardiomyocytes (ten Tusscher - TNNP) and a database of human physiological, genotypic and demographic data enabling the prediction of the QT prolongation in humans based on the in vitro data [7, 22]. To account for the heterogeneities in ionic currents between endocardial, midmyocardial and epicardial cells 1D fibre paced at the epicardial side was constructed. The 50:30:20 distribution of the endo-, mid- and epicardial cells was used together with a diffusion coefficient equal to 0.0016 cm2/ms. The Forward Euler method was used to integrate model equations. Integration results were used to calculate a pseudo-ECG. First and last QRS were excluded from the pseudo-ECG. A space step and a time step were set to Δx=0.01 mm and Δt=0.01 ms, total simulation time was set to 10,000 ms.
To account for the drugs triggered ionic currents modifications a specific equation describing the current of interest was multiplied by the inhibition factor accordingly to the in vitro values provided by the literature search, which described the concentration dependent ionic current inhibition. The inhibition factor was calculated with the use of the Hill equation [Equation 1].
IC50 - concentration responsible for the 50% inhibition of the ionic current
n - Hill equation parameter
DRUG CONCENTRATION - active drug concentration [μM]
The population variability of other parameters was mimicked by applying the virtual population generator as described previously [23, 24]. The circadian heart rate variability was introduced into a simulation by the use of the model described below.
3.8 μM of MOXI was used for the operational concentration, mimicking the average maximum free plasma concentration (Cmax) after a 400 mg oral dose by correcting the total Cmax concentration obtained from the available literature with human plasma protein binding [25, 26]. The IC50 values for various cardiac ion channels were obtained from the tox-database.net system and are presented in Table 1 . For the simulation studies only IKr data was used (with Hill equation parameter n = 1) since for the tested MOXI concentration INa and ICa inhibition were negligible.
The ToxComp simulations, at a MOXI concentration equal to Cmax at different times of the day (4.00, 8.00, 12.00, 16.00, 20.00, and 24.00), were performed in triplicate trials of 20 individuals (total 3×20=60) with an equal number of male and female subjects (10/10 in each study). Baseline QT and QTcF (QT interval corrected by heart rate using Fridericia correction method ) were also obtained for each subject at the above defined time points of the day by simulating a response at MOXI concentration equal to zero. This allowed us to apply two baseline correction formulas i.e., (1) Single point baseline correction when the average baseline QTcF estimates from all subjects was used to obtain baseline corrected QTcF (ΔQTcF) and (2) Individualised baseline correction (ΔQTcFi) when baseline QTcF at a given time of day for each individual subject was used to calculate ΔQTcF. The developed circadian effect model was validated by comparing the simulated results with respective clinical outcomes obtained using single point and individualised baseline correction formula.
The aim of this section is to describe the steps undertaken to create a multivariate linear regression model of the relationship between RR (dependent variable) and a set of independent variables which were comprised of Age, Sex and Time of the day (abbreviated as Hour). In addition to the model development, the simulation of RR values from the model was given attention. During the model development phase it was found that the distribution of RR data was mildly skewed to the right. This issue was dealt with the use of log-transformation of RR measurements.
The preliminary model included dummy variables for each Hour, dummy variable for Sex, quadratic effect of Age and all pairwise interactions. Many terms in this preliminary model were found not significant, as indicated by the robust t-tests for model coefficients. The robust t-test takes into account the dependence between observations coming from the same subject. Non-significant terms (p-value>0.05) were dropped from the model sequentially one by one. Additionally it was found that the Hour variable may be modelled more parsimoniously by a linear combination of sine and cosine functions. The final model received the following form:
β – regression coefficient,
ε – normally distributed error term,
σ – standard deviation of error term,
Sex – 1 for males, 0 for females,
Age – age in years,
Hour – value from 0–24 range.
The model parameters were estimated using the maximum likelihood method in the R system for statistical computing and are given below :
All coefficients were statistically significant. Figure 1 shows the graphical representation of the model.
The visual predictive checks showed an overall good fit of the model to the data. The dispersion of residuals did not exhibit dependence upon any of the explanatory variables. The coefficient of determination R2 = 0.39, i.e. 39% of the variation observed in log(RR) can be explained by the estimated regression equation.
For the PhysioNet data set RMSE (Root Mean Squared Error) = 128 ms and MAPE (Mean Absolute Percentage Error) = 12.3%. In the case of the validation data set RMSE = 165 ms and MAPE = 17.1%. Figure 2 shows the representative two best and two worst fitted cases from the validation data set.
It is a well-known fact that estimators derived by the method of maximum likelihood have many desirable properties, on the other hand ML estimators are sensitive to the presence of outliers in the data. In order to check the stability of estimates the regression parameters were additionally estimated using a wide range of, so called, robust regression estimators. Parameter estimates were virtually unchanged from the maximum likelihood ones.The residuals (
) within each subject showed strong temporal dependence. Consequently the stationary autoregressive model of the following form was fitted:
α – autoregression coefficient
η – normally distributed error term
τ – standard deviation of error term
P – order of autoregression
It was found that only a large value of P (P = 180) was able to account for the long memory observed in the data and gives a satisfactory fit. The estimate of the standard deviation of error term was 0.096.
The following algorithm facilitates the generation of random RR values from the estimated model:
Generate random sequence from equation [Eq. 4].
Choose the values of explanatory variables and using [Eq. 3] compute at the selected time points.
Add the computed to values generated in Step 1 at the selected time points, exponentiate the result.
An electronic supplement (Additional file 1) to the article contains a fully functional Excel implementation of the model.
Model usability results
QT response at a MOXI concentration of 3.8 μM at different times of the day was simulated using the ToxComp platform. When single point baseline correction was applied, the ΔQTcF obtained at different times of the day at a MOXI concentration of 3.8 μM is shown in Figure 3. Generally, during TQT studies an oral dose of MOXI 400 mg is given in the morning, hence Cmax is obtained between morning and noon i.e. between 8.00 to 12.00. At 8.00 and 12.00, the simulated ΔQTcF were respectively 10.75 ms and 15.87 ms, which is consistent with the clinically observed ΔQTcF of 10–15 ms . When an individualized baseline with circadian correction was applied, ΔQTcFi of 5–6 ms was obtained at all times of day which is also consistent with the clinically observed ΔQTcFi of 6–7 ms .
The obtained results indicate that the implemented heart rate circadian rhythm model within the ToxComp platform is able to simulate the circadian variability of a model drug. Comparison with the published clinical trial results for the model drug suggest that the ToxComp system, with the built-in heart rate model, is able to realistically represent cardiac electrophysiological drug effect with its variability, regardless of the correction method. However, as there are a large number of parameters describing the physiological variability in the virtual population generator ToxComp module, it is hard to precisely assess the net influence of the heart rate on the system output, and a separate study will be run for such a purpose. Although results of the simulations run without the use of a heart rate model, i.e. with the use of a constant value of HR during the day, suggested that an important component is missing (internal not-published results).
The goodness-of-fit metrics indicated a good fit to the data, e.g. Mean Absolute Percentage Error equals to 12.3%. It could, although, suggest that the model overfits the data. This was however contradicted by the validation results with MAPE equals to 17.1%. The result was obtained despite the differences in the age and sex distribution of the subjects between the training and validation data sets.
The main limitation of the model is the size of the training data set. Despite the fact that every subject gave rise to thousands of measurements, the number of independent subjects is still only 18. The future research will aim to combine data from experiments which acquire sparse RR data, e.g. every 20–40 minutes, with dense PhysioBank data.
It can be concluded that a flexible and practically useful model describing the heart rate circadian variation has been developed. Its structure allows for easy implementation, which is greatly facilitated by the provided electronic supplement, and its use in simulation studies.
Chinese Hamster Ovary cells
Human Embryonic Kidney cells
In vitro–in vivo extrapolation
Mean Absolute Percentage Error
QT interval corrected for heart rate with the Fridericia correction
Difference between the average of baseline QTcF of all subjects and QTcF in the presence of a drug
Difference between the baseline QTcF of each subject at a given time of day and QTcF in the presence of a drug for the same subject
Coefficient of determination
Root Mean Squared Error
Interval from the onset of one QRS complex to the onset of the next QRS complex
Torsade de Pointes
Thorough QT study
Hanton G: Preclinical cardiac safety assessment of drugs. Drugs R&D. 2007, 8 (4): 213-228.
Inoue H, Yamanaka S: The use of induced pluripotent stem cells in drug development. Clin Pharmacol Ther. 2011, 89: 655-661. 10.1038/clpt.2011.38.
ICH 2005 International Conference on Harmonization of Technical Requirements for Registration of Pharmaceuticals for Human Use: The Clinical Evaluation of QT/QTc Interval Prolongation and Proarrhythmic Potential for Non-Antiarrhythmic Drugs (ICH E14). 2005, Rockville, Md: US Department of Health and Human Services Food and Drug Administration
Redfern WS, Carlsson L, Davis AS, Lynch WG, MacKenzie I, Palethorpe S, Siegl PK, Strang I, Sullivan AT, Wallis R, Camm AJ, Hammond TG: Relationships between preclinical cardiac electrophysiology, clinical QT interval prolongation and torsade de pointes for a broad range of drugs: evidence for a provisional safety margin in drug development. Cardiovasc Res. 2003, 58 (1): 32-45. 10.1016/S0008-6363(02)00846-5.
Yao X, Anderson DL, Ross SA, Lang DG, Desai BZ, Cooper DC, Wheelan P, McIntyre MS, Bergquist ML, MacKenzie KI, Becherer JD, Hashim MA: Predicting QT prolongation in humans during early drug development using hERG inhibition and an anaesthetized guinea-pig model. Br J Pharmacol. 2008, 154 (7): 1446-1456. 10.1038/bjp.2008.267.
ToxComp.net. www.tox-comp.net accessed 07-02-2013
Polak S, Wisniowska B, Fijorek K, Glinka A, Polak M, Mendyk A: ToxComp ‐ in vitro–in vivo extrapolation system for the drugs proarrhythmic potency assessment. Com Cardiol. 2012, 39: 789-793.
Polak S, Wisniowska B, Fijorek K, Glinka A, Polak M, Mendyk A: Modeling and simulation approach for assessing proarrhythmic potency of the non‐cardiological drugs. Com Cardiol. 2012, 39: 761-765.
Polak S, Wisniowska B, Fijorek K, Glinka A, Polak M, Mendyk A: In silico prediction of the drug overdose consequences at the heart electrophysiology level. Com Cardiol. 2012, 39: 889-893.
Aschoff J: Circadian Rhythms in Man. A self-sustained oscillator with an inherent frequency underlies human 24-hour periodicity. Science. 1965, 148 (3676): 1427-1432. 10.1126/science.148.3676.1427.
Kanabrocki EL, Sothern RB, Scheving LE, Vesely DL, Tsai TH, Shelstad J, Cournoyer C, Greco J, Mermall H, Ferlin H: Reference values for circadian rhythms of 98 variables in clinically healthy men in the fifth decade of life. Chronobiol Int. 1990, 7 (5–6): 445-461.
Piotrovsky V: Pharmacokinetic-pharmacodynamic modeling in the data analysis and interpretation of drug-induced QT/QTc prolongation. AAPS J. 2005, 7 (3): E609-E624. 10.1208/aapsj070363.
Nakagawa M, Iwao T, Ishida S, Yonemochi H, Fujino T, Saikawa T, Ito M: Circadian rhythm of the signal averaged electrocardiogram and its relation to heart rate variability in healthy subjects. Heart. 1998, 79 (5): 493-496.
Massin MM, Maeyns K, Withofs N, Ravet F, Gerard P: Circadian rhythm of heart rate and heart rate variability. Arch Dis Child. 2000, 83: 179-182. 10.1136/adc.83.2.179.
Li X, Shaffer ML, Rodriguez-Colon S, He F, Wolbrette DL, Alagona P, Wu C, Liao D: The circadian pattern of cardiac autonomic modulation in a middle-aged population. Clin Auton Res. 2011, 21 (3): 143-150. 10.1007/s10286-010-0112-4.
Bonnemeier H, Richardt G, Potratz J, Wiegand UK, Brandes A, Kluge N, Katus HA: Circadian Profile of Cardiac Autonomic Nervous Modulation in Healthy Subjects: Differing Effects of Aging and Gender on Heart Rate Variability. J Cardiovasc Electrophysiol. 2003, 14: 791-799. 10.1046/j.1540-8167.2003.03078.x.
Eisenbruch S, Harnish MJ, Orr WC: Heart rate variability during waking and sleep in healthy males and females. Sleep. 1999, 22 (8): 1067-1071.
Cavallari JM, Fang SC, Mittleman MA, Christiani DC: Circadian variation of heart rate variability among welders. Occup Environ Med. 2010, 67 (10): 717-719. 10.1136/oem.2010.055210.
Cugini P, Bernardini F, Cardarello CM, Coda S, Curione M, De Francesco GP, De Rosa R, Dutto L, Fontana S, Mammarella A, Paoletti V, Paradiso M, Pellegrino AM: Circadian rhythm of heart rate in myotonic dystrophy. J Clin Basic Cardiol. 2000, 3 (3): 181-186.
D'Negri CE, Marelich L, Vigo D, Acunzo RS, Girotti LA, Cardinali DP, Siri LN: Circadian periodicity of heart rate variability in hospitalized angor patients. Clin Auton Res. 2005, 15 (3): 223-232. 10.1007/s10286-005-0280-9.
Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PC, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE: PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals. Circulation. 2000, 101 (23): 215-220. 10.1161/01.CIR.101.23.e215.
ten Tusscher KH, Noble D, Noble PJ, Panfilov AV: A model for human ventricular tissue. Am J Physiol Heart Circ Physiol. 2004, 286 (4): H1573-H1589.
Polak S, Fijorek K, Glinka A, Wisniowska B, Mendyk A: Virtual population generator for human cardiomyocytes parameters. In silico drug cardiotoxicity assessment. Toxicol Mech Methods. 2012, 22 (1): 31-40. 10.3109/15376516.2011.585477.
Polak S, Fijorek K: Inter-individual variability in the preclinical drug cardiotoxic safety assessment - analysis of the age - cardiomyocytes electric capacitance dependence. J Cardiovasc Transl Res. 2012, 5 (3): 321-332. 10.1007/s12265-012-9357-8.
Grosjean P, Urien S: Reevaluation of moxifloxacin in pharmacokinetics and their direct effect on the QT interval. J Clin Pharmacol. 2012, 52: 329-338. 10.1177/0091270011398361.
Stass H, Dalhoff A, Kubitza D, Schuhly U: Pharmacokinetics, safety, and tolerability of ascending single doses of Moxifloxacin, a new 8-methoxy Quinolone, administered to healthy subjects. Antimicrob Agents Chemother. 1998, 42 (8): 2060-2065.
Polak S, Wisniowska B, Glinka A, Polak M: Tox-database.net: a curated resource for data describing chemical triggered in vitro cardiac ion channels inhibition. BMC Pharmacol Toxicol. 2012, 13 (6): 1-13.
Alexandrou AJ, Duncan RS, Sullivan A, Hancox JC, Leishman DJ, Witchel HJ, Leaney JL: Mechanism of hERG K+ channel blockade by the fluoroquinolone antibiotic moxifloxacin. Br J Pharmacol. 2006, 147 (8): 905-916. 10.1038/sj.bjp.0706678.
Champeroux P, Ouillé A, Martel E, Fowler JS, Maurin A, Richard S, Le Guennec JY: A step towards characterisation of electrophysiological profile of torsadogenic drugs. J Pharmacol Toxicol Methods. 2011, 63 (3): 269-278. 10.1016/j.vascn.2011.01.001.
Fridericia LS: Die sytolendauer in elektrokardiogramm bei normalen menschen und beiherzkranken. Acta Med Scand. 1920, 53: 469-486.
R Development Core Team: R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing. 2012, Vienna, Austria: , http://www.R-project.org/ , 3-900051-07-0,
Garnett C, Beasley N, Bhattaram VA, Jadhav PR, Madabushi R, Stockbridge N, Tornoe CW, Wang Y, Zhu H, Gobburu JV: Concentration-QT relationships play a key role in the evaluation of proarrhythmic risk during regulatory review. J Clin Pharmacol. 2008, 48: 13-18. 10.1177/0091270007307881.
Project financed by Polish National Center for Research and Development. LIDER project number LIDER/02/187/L-1/09. KF acknowledges financial support from The Foundation for Polish Science.
The authors declare that they have no competing interests.
KF and NP participated equally in the study conductance, KF – model development, NP – model validation. Study idea and methodology development come from SP and KF. KF, NP and SP were equally involved in the manuscript writing. ŁK, KS-S, KK-J were responsible for the clinical data gathering and control from the clinical point of view. All authors read and approved the final manuscript.
Electronic supplementary material
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.