A myofibre model for the study of uterine excitation-contraction dynamics

As the uterus remodels in preparation for delivery, the excitability and contractility of the uterine smooth muscle layer, the myometrium, increase drastically. But when remodelling proceeds abnormally it can contribute to preterm birth, slow progress of labour, and failure to initiate labour. Remodelling increases intercellular coupling and cellular excitability, which are the main targets of pharmaceutical treatments for uterine contraction disorders. However, the way in which electrical propagation and force development depend on intercellular coupling and cellular excitability is not fully understood. Using a computational myofibre model we study the dependency of electrical propagation and force development on intercellular coupling and cellular excitability. This model reveals that intercellular coupling determines the conduction velocity. Moreover, our model shows that intercellular coupling alone does not regulate force development. Further, cellular excitability controls whether conduction across the cells is blocked. Lastly, our model describes how cellular excitability regulates force development. Our results bridge cellular factors, targeted by drugs to regulate uterine contractions, and tissue level electromechanical properties, which are responsible for delivery. They are a step forward towards understanding uterine excitation-contraction dynamics and developing safer and more efficient pharmaceutical treatments for uterine contraction disorders.

isolated cells or simplified depolarization and contraction waves at the whole tissue level. The cellular level models [12][13][14] are useful in studying AP morphologies and calcium dynamics. However, these models exclude intercellular dynamics and therefore cannot be used to study electrical conduction and consequent mechanical activity.
Whole organ models, on the other hand, are useful for studying geometrical trends of excitation and contraction, and for solving the forward and inverse problems of uterine electrophysiology [15][16][17][18] . These models, however, require drastic reductions of the cellular models and can miss dynamically regulated cellular responses.
In this study, we developed a myofibre model to study the regulatory mechanisms of electrical conduction, conduction block, and force development. This study characterizes the propagation of the electromechanical wave in the myometrium. Our results introduce a theoretical framework that is useful to interpret experimental observations. Moreover, our model bridges between existing cellular level models and whole organ models. This bridge is necessary to understand the interplay between the cellular targets of drugs developed to regulate uterine contractions, and the tissue and organ level electromechanical properties responsible for delivery.
Our model consists of a series of cells, described by a detailed ionic current model, coupled longitudinally. The intercellular interactions are described by a system of equations that governs the electromechanical wave propagation.
We simulated the propagation of an electrical stimulus along the myofibre to study the effects of intercellular coupling and cellular excitability on uterine electrophysiology and force development. We made four observations: (1) Intercellular coupling determines the conduction velocity (CV) and remodels the AP. (2) Force development has a binary dependence on intercellular coupling: the force developed is constant under high intercellular coupling, and only a passive tension is developed under low intercellular coupling. (3) Increased cellular excitability can recover conduction that was blocked by low intercellular coupling, but it cannot recover the reduction in CV. Lastly, (4) Force development has a discontinuous dependency on cellular excitability.
Our results contribute a step forward in understanding the responses to tocolytic and uterotonic drugs. The mechanisms we explored, remodel in preparation for parturition. These mechanisms, in turn, are targeted by all the available uterotonic and tocolytic drugs. Consequently, better understanding of these mechanisms is necessary to develop more efficient and safer drugs.

Results
Model overview. We model a myofibre as a longitudinal concatenation of cells. Our multiscale model combines a detailed USMC model with a novel electromechanical myofibre model (Fig. 1a). Each cell in the myofibre includes a detailed ionic current model, calcium handling dynamics, actin-myosin interactions, and a forceproducing mechanism. The myofibre model consists of 70 cells serially connected. We chose to model 70 cells in the myofibre since this length is sufficient to contain the entire depolarization and contraction waveforms, yet not too long to make the simulations computationally expensive. To exclude end effects, the first and last 10 cells of the myofibre are not considered for analysis.
The intercellular interactions are modeled using a system of differential equations, the electrical propagation is modeled through the cable equation, while the mechanical transduction is modeled by a system of 3n equations, where n is the number of cells in the myofibre. The electrical activity modulates the mechanical activity through the calcium dynamics, actin-myosin cycling, and cellular contraction model. Conversely, the mechanical activity modulates the electrical activity by varying the geometrical parameters in the cable equation. In particular, the mechanical wave arises from contracting and relaxing cells through the myofibre, which have varying lengths depending on their contractile state. This effect is taken into account through the spatial component of the wave equation. The complete formulation of our model appears in the "Methods" section.
We stimulate the myofibre from one end with a 30 ms long depolarization stimulus of − 5 pA/pF. This stimulus produces a single AP that propagates along the myofibre. Although APs in the contracting myometrium usually appear in bursts, the propagation of a single AP is useful to study the mechanisms responsible for the excitability and contractility of the myometrium. Considering a single AP simplifies the analysis by excluding the effects related to the properties of the burst, such as the firing rate. Additionally, many parameters used to quantify excitability and contractility, such as conduction velocity and action potential duration, are not well defined for a burst of APs.
Our model recapitulates experimental observations. The single cell model is based on the models by Tong et al. 12 and Testrow et al. 13 , as well as on the previous models by Yang et al. 19 and Hai and Murphy 20 . The developed myofibre model is based on Newton's laws and recapitulates physiological behaviours. The AP waveforms obtained have physiological morphologies 21 , and propagate with experimentally observed conduction velocities 22 (Fig. 1b). Additionally, the intracellular calcium transients (Fig. 1c), cellular contraction dynamics (Fig. 1d), and contractile force development (Fig. 1e) behave as expected, based on previously reported experimental recordings 12,20,21,23,24 . CV and AP morphology dependence on intercellular coupling. Intercellular coupling is achieved through gap junctions (GJ), which form channels that connect the cytosol of two adjacent cells and allow the flow of ions and small molecules. GJs consist of two connexons, each formed by six connexins, anchored to each of the connected cell's membranes. Connexin 26,32,37, and 43 form myometrial GJs 25 . The expression level of GJs, as well as the type of connexins used to form those GJs, regulate the intercellular resistivity (IR) 26 .
The healthy myometrium remodels to induce parturition 25,27 . During pregnancy, connexins are expressed at a basal level, providing minimal intercellular coupling and preventing electrical conduction. Before birth, connexin expression and GJ formation is upregulated. Intercellular coupling increases, allowing APs to propagate across Scientific RepoRtS | (2020) 10:16221 | https://doi.org/10.1038/s41598-020-72562-x www.nature.com/scientificreports/ cells. This process, however, may go amiss, contributing to either premature contractions or slowly progressing labours 28 . We simulated electrical conduction along a myofibre under different intercellular coupling conditions. In our model, intercellular coupling is determined by the IR between adjacent cells, and intercellular coupling is inversely related to IR. We varied the IR between the cells in the myofibre (as detailed in the "Methods" section) and determined its effects on conduction velocity, cellular depolarization, and conduction block.
Conduction velocity is defined as the propagation velocity of the depolarizing wave through the myometrium. This velocity has been reported to usually range between 1 and 3 cm/s in-vivo, but it has also been reported to be as high as 12 cm/s 22 . Cellular depolarization is driven mostly by a rapid calcium influx and is characterized by its maximal upstroke slope max t ∂v m ∂t and the duration of its upstroke. These two parameters define the ability of a USMC to propagate an AP. In our simulations, the time difference between the depolarization of adjacent cells is much shorter than the upstroke duration, so the upstroke section of the AP drives the propagation. The duration of the upstroke is measured from when the cell achieves 5% of the maximal depolarization until it achieves its maximal transmembrane voltage. Conduction block is defined as the inability of the stimulus to traverse the entirety of the myofibre.    www.nature.com/scientificreports/ Electrical conduction behaves differently for high and low intercellular coupling, modeled by IR values above and below 25 cm, respectively. In the high intercellular coupling domain, CV is very sensitive to IR value changes, whereas CV has a low sensitivity for IR in the low intercellular coupling domain (Fig. 2a). Conduction block first occurs at 200 cm. For IR values above that threshold, a stimulus propagates to only a few cells. From the reported data on CV and our simulations, we infer that the IR in excitable myometrium ranges between 10 and 100 cm.
As intercellular coupling increases, membrane depolarization becomes faster. The maximal upstroke slope increases with increasing intercellular coupling (Fig. 2b), while the upstroke duration decreases with increasing intercellular coupling (Fig. 2c). As with conduction velocity, the sensitivity of the upstroke slope, and its duration, to intercellular coupling is high in the high intercellular coupling domain and low in the low intercellular coupling domain.
Force development dependence on intercellular coupling. As long as conduction is unblocked, the tensile force developed by the myofibre model is insensitive to IR. At every time step, our model calculates the tensile force developed by the myofibre in response to an initial stimulus. The forces developed under different IR values are compared by their amplitudes and normalized cumulative force values, defined as the integral of the force over 14 s after the initial stimulus, normalized by the passive force (i.e., without stimulation) developed by the myofibre over the same time period. This time period is the shortest time for which all the simulations return to their initial state after stimulation. We calculated the normalized cumulative force (Fig. 3a), as well as the amplitude of the force developed (Fig. 3b), for different IR values. The force developed has a weak dependence on IR for low IR values, but drops sharply around 200 cm , where conduction block occurs. For higher IR values, the normalized area under the curve approaches unity; there is no force developed besides passive tension.
The robustness of the force developed against IR arises from ionic channel dynamics. The AP adapts to changing conduction velocities, buffering the calcium influx and consequently conferring insensitivity to the force developed from IR. The cellular dynamics of a cell in the middle of the myofibre explain this phenomenon. We consider the AP, the calcium influx, and the force-producing mechanism for the 30th cell in a 70-cell-long myofibre. The myofibre is stimulated as mentioned in the model overview, and the IR is varied between 10 and 400 cm. For IR below 200 cm, the entire myofibre is excited, at 300 cm the cell selected is in the area just before conduction block occurs, and at 400 cm the cell is past the conduction block.
As the IR increases, conduction velocity decreases, and therefore the AP at the 30th cell is progressively delayed. The waveform amplitude is largely decreased just before conduction block, and the cell fails to activate for 400 cm (Fig. 4a). The waveforms obtained for normal propagation are compared by the action potential duration at 90% repolarization ( APD 90 ), defined as the time required for the cell to restore to 90% its membrane potential, and through the AP amplitude. As the IR increases, the APD 90 increases while the amplitude decreases (Fig. 4b). These two trends work in tandem to buffer the force driving calcium inward.
The intracellular calcium trace is characterized by a sharp upstroke followed by a slow inactivation (Fig. 4c). The increased IR delays the onset of the upstroke, in agreement with the AP delay. The intracellular calcium trace pattern responds to the AP dynamics: as the IR increases, the calcium influx current's amplitude decreases, while it's duration prolongs (Supplementary Fig. S2). These balancing forces keep the decay slope unchanged    (Fig. 4d). By controlling the cross-bridge cycling dynamics, the intracellular calcium concentration regulates the force developed. Detailed descriptions of this mechanism can be found in the review by Aguilar et al. 29 . Briefly, myosins cycle through multiple states, and as they cycle, they attach and detach from actin filaments, producing a contractile force on the cell. The fraction of myosin bound to actin determines the force developed (Fig. 4e). In smooth muscles, the cycling of myosin through its multiple states is regulated by the myosin light chain kinase (MLCK), and MLCK's activity, in turn, is regulated by the intracellular calcium concentration. This system was modeled by Yang et al. 19 and included in our model (as detailed in the "Methods" section). The fraction of myosin in the force-producing states over time follows the intracellular calcium dynamics.
The force developed by the myofibre is calculated through our system of coupled equations (see "Model overview" and "Methods"). The bound myosin fractions of the cells in the myofibre determine the force developed by the myofibre (Fig. 4f). With increasing IR values, the force upstroke losses its sharpness and the peak force is attained at a later time. The buffered cross-bridge cycling dynamics maintain the force amplitude and normalized cumulative force for IR values below 150 cm. This buffer effect is lost for IR values above 150 cm, which is expected since conduction block occurs at 200 cm (Fig. 3).
Electrical conduction dependence on cellular excitability. Uterine smooth muscle cell excitability is mainly supported by calcium currents. The AP initial depolarization is driven by a long-lasting calcium current (I CaL ), with a secondary transient calcium current (I CaT ) and sodium current (I Na ) 12,30,31 . This configuration differs from the well-studied cardiac and neural systems, where a single rapid sodium current is responsible for depolarizing the cell. Calcium channels in USMC require a sustained stimulus to activate and initiate the AP. The depolarizing stimulus for a given cell is provided by its adjacent depolarized cell through the GJ connecting both cells.
USMC excitability is highly regulated, but the effect of cellular excitability modulation on electrical conduction is incompletely understood. Ion-channels responsible for USMC excitability are upregulated as parturition approaches. Common uterotonic drugs, such as oxytocin, enhance cellular excitability to induce slowly progressing labours. On the other hand, tocolytic drugs are ion-channel blockers that aim to inhibit cellular excitability 32 .
To study the dependency of electrical conduction and force development on intercellular coupling and cellular excitability, we simulated the AP propagation through a myofibre with varying IR values, while modulating the cell's conductivity for either I CaL , I CaT , or I Na (Fig. 5).  www.nature.com/scientificreports/ Conduction block caused by low intercellular coupling is overcome with increasing cellular excitability. Increasing the conductance of the ion channels responsible for the I CaL (Fig. 5a), I CaT (Fig. 5b), or I Na (Fig. 5c) currents restores conduction for IR values above the threshold for conduction block, and the cellular excitability upregulation required to restore conduction increases with decreasing intercellular coupling. Conduction block is most sensitive to I CaL currents: small I CaL downregulations block conduction even in well coupled myofibres, while small upregulations overcome blockage in poorly coupled myofibres.
Although conduction blockage is restored with increased cellular excitability, the conduction velocity is not. With increased cellular excitability, the conduction velocity of poorly coupled myofibres is considerably smaller than that of well coupled myofibres with normal cellular excitability. For example, increasing the conductivity of I CaL in myofibres with IR above 100 cm leads to only marginal increases in the conduction velocity (Fig. 5a). Within the physiological range of our model, CV is less sensitive to cellular excitability changes than to intercellular coupling changes.
Force development dependence on cellular excitability. The effect of cellular excitability modulation on force development depends on the type of electrical current. In this study we explored the dependency on the main depolarizing currents, i.e., I CaL , I CaT , and I Na . As the conductivity for I CaL increases, but does not overcome the conduction block, the force developed remains at the passive tension level. Then, as the conduction block is overcome, the force developed increases rapidly. After that, the force continues to increase gradually with increasing I CaL conductivity (Fig. 5d).
The force developed increases for the entire range of I CaT conductivity modulation (Fig. 5e). As conduction block is overcome, the plot of the force developed with respect to I CaT modulation switches from the lower to the upper curve.
The force developed has a binary behavior with respect to I Na conductivity modulation (Fig. 5f). The force remains at the passive tension level while conduction is blocked, then takes a higher value after the conduction block is overcome. Further increasing the I Na conductivity, however, does not lead to increasing tensile forces.
For all current types, the force developed before and after overcoming the conduction block is independent of the intercellular coupling. Lower intercellular coupling implies that a higher cellular excitability increase is required to overcome the conduction block. However, once the conduction block is overcome, the force developed depends on I CaL and I CaT conductivity, and not on intercellular coupling.

Discussion
We developed a uterine myofibre model and used it to explore electrical conduction and force development. This model integrates validated models of membrane currents, calcium dynamics, and force development in a novel mathematical framework. A coupled system of differential equations governs the propagation of the electromechanical wave along a series of cells modelling a myofibre. To study the regulation of electrical conduction and Our model simulation showed that intercellular coupling is the main regulator of CV. Very low intercellular coupling induces conduction block, and for higher intercellular coupling we identified two domains. In the low intercellular coupling domain, we observed a slow CV that is insensitive to intercellular coupling changes. On the other hand, in the high intercellular coupling domain, we observed a high CV that is also very sensitive to intercellular coupling. Although cellular excitability changes, simulated by modulation of the ion channels' conductance, can recover conduction block caused by low intercellular coupling, the effect of those changes on CV is marginal. Our simulations also showed that the force developed is buffered for varying intercellular coupling values, provided the stimulus propagates along the entire myofibre. The simulated traces for the voltage, calcium, and myosin state explain this compensatory mechanism. Lastly, we explored the effect of cellular excitability on force development and found that upregulating the conductivity of I Na and I CaL has a minimal effect on the force developed, once conduction is established. Upregulating I CaT , on the other hand, has a profound effect on the force developed.
Our results contribute to a better understanding of uterine electrophysiology. To the best of our knowledge, this is the first computational study to explore and characterize the propagation of the electromechanical wave in the myometrium. We studied the effects of intercellular coupling and cellular excitability on electrical propagation and force development. Intercellular coupling and cellular excitability are remodeled at the end of pregnancy, in preparation for parturition. Moreover, all the available tocolytic and uterotonic drugs target these mechanisms. Therefore, the exploration of these mechanisms is fundamental to understanding the remodelling processes of the myometrium and to developing more efficient and safer drugs to modulate them.
All models, including this one, have limitations. Although multiple AP morphologies have been reported in the myometrium 33 , we considered only a single AP morphology. Our simulations are based on the spiked AP elicited by a brief depolarization of a USMC. Tong et al. 12 showed that their USMC model, which was used in this study, can reproduce multiple AP morphologies and the effects of AP trains. Our study explores the intrinsic properties of the myofibre but not the dependency on the AP morphology. Additionally, multiple ionic current models of USMC are available. We chose to use the model of Tong et al. based on their experimental validation and a preliminary stability analysis. Our framework, however, is modular, and the USMC ionic current model can be replaced with newer models.
An important limitation of our model is that it has been developed from available smooth muscle models, rather than from dedicated experimental observations. Firstly, the models integrated in our study were developed from experimental observations recorded by several groups using different experimental systems and species. Secondly, the models used to describe myosin cycling 20 and the subsequent cellular contraction 19 , were developed for vascular smooth muscles. In the future, as experimental and theoretical uterine electrophysiology advances, uterine excitation-contraction models will preferably be developed based on consistent experimental observations and theoretical frameworks.
Future extensions of our study include investigating other regulatory mechanisms of conduction and force development in the myometrium. In particular, we note that the myometrium expresses stretch activated potassium channels. These channels open in response to stretch and inhibit membrane excitability. Multiple groups have shown that these ion channels are expressed during pregnancy and downregulated towards term. This suggests that these ion channels play a role in maintaining the quiescent electrical state of the uterus during pregnancy 11,34,35 .
Moreover, other mechanosensitive ion channels have been found in the myometrium, e.g., swell-activated chloride channels, non-specific cation channels, and calcium channels 11 . Our model implements the ion channel model developed by Tong et al. 12 but does not incorporate mechanosensitive ion channels. In the future, this ion channel model could be replaced for another that includes mechanosensitive ion channels to study in the electromechanical feedback in further detail.
Additionally, we are interested in the effects of myosin light chain kinase regulation. This kinase, regulated by oxytocin and progesterone, is thought to be a key regulator of myometrial excitability and contraction 10,36 . Lastly, while uterine APs usually appear in bursts and our model can simulate the bursting type AP 12,13 , as shown in Supplementary Fig. S3, we only used single spiked APs for our simulations. Using single APs simplifies the analysis of the mechanisms underlying the excitability and contractility of the myofibre. However, additional complex observations can be made from analyzing the propagation of the bursting type AP. Further study of these mechanisms may help prevent preterm birth and assist in labor induction.

Methods
USMC model. We model a myofibre as a longitudinal concatenation of USMC ionic models. Each cell is represented by the detailed ionic current model developed by Tong et al. 12 . The complete formulation and implementation of this model can be found in the original publication. Briefly, the ionic currents, the transmembrane voltage, and the intracellular calcium concentration are related through a system of over 100 differential equations. This model contains 15 different ionic currents, most of which are modeled through the Hodgkin Huxley formalism, i.e.,

Scientific RepoRtS
| (2020) 10:16221 | https://doi.org/10.1038/s41598-020-72562-x www.nature.com/scientificreports/ where I is the transmembrane ion current, ḡ is the maximal conductance for the given current, y i are the gating variables, v is the membrane potential, E rev is the Nernst potential for the ion X, R is the gas constant, T is the absolute temperature, F is the faraday constant, [X] o and [X] i are the extracellular and intracellular ion concentrations, and y ∞ i and τ y i are the steady state and time constants for channel i and depend on the transmembrane voltage and intracellular calcium concentration.
Excitation model. A depolarizing current is injected into the first cell in a series and triggers an AP. The depolarization of USMCs results mostly from inward calcium currents and promotes the depolarization of the following cell through their electrical coupling. The resulting depolarization wave is governed by the cable equation. Assuming an homogenous extracellular medium, the cable equation is given by 37 where V(x, t) is the transmembrane voltage along the cable; I ion (x, t) is the ionic current based on the model of Tong et al.; I stimulus (x, t) is the stimulatory current, which in our case is spatially constrained to the first cell; C m is the specific membrane capacitance; σ is the intercellular conductivity; and χ is the cell's surface to volume ratio. We assume that the cells' volume and surface are deformable but incompressible. In other words, while the cells' shape can change, the total volume and surface area remain constant. Therefore, in our model the surface to volume ratio (χ) remains constant.
The intercellular coupling is the product of the intercellular cross sectional area, A ics , that is the cross sectional area at the intercellular junction, and the intercellular conductivity, σ . The intercellular coupling is inversely proportional to IR, denoted by ρ.
In our model, we assume that the intercellular cross sectional area, A ics , remains constant throughout the simulation. This assumption is supported by the fact that the intercellular junctions are rich in anchoring protein complexes that hold the membrane's surface area at the intercellular junction. We do not assume, however, that the cellular cross sectional area remains constant throughout the cell. On the contrary, the cellular cross sectional area can change, under the constraints that the total volume, total surface area, and cross sectional area at the intercellular junction remain constant. Here, nm is the Hill coefficient; CaMLCK is the half-saturation concentration of MLCK; and K 2 , K 3 , K 4 , K 5 , and K 7 are myosin state transition constants.

Cellular contraction model. [Ca
We modeled the contractile mechanics of a single cell using the model developed by Yang et al. 19 . In this model, the force is developed by a passive element connected in parallel to three serially connected force-generating elements. The force-generating elements are the cross bridge elasticity, the active force-generating element, and the series viscoelastic element. A listing of the mechanical parameters appears in Table 1.
The passive force (p) is given by The force generated by the cross bridge elasticity (x) is given by The tensile force generated by the cell is given by Myofibre contraction model. For a travelling AP in a myofibre of i = 1 . . . n cells, the myosin cycling is unique to each cell and is dependent on its electrical activation. Our problem, then, is as follows: given AM i (t) and AMp i (t) , we are interested in finding the lengths of each cell in the myofibre and the tensile force developed by the entire myofibre. To solve this problem, we need to find the lengths of the cellular compartments of each cell at every time step. For a myofibre of n cells, we get 3n unknowns: l i j (t) i=1...n j=a,s,x . Combining Eqs. (9), (10), and (13) we get Combining Eqs. (9), (11), and (13), we then get (12) l c = l a + l x + l s .  (14) and (17), we get The last equation needed to solve the system comes from the isometric condition. The myofibre's total length, but not an individual cell's length, is held constant. This condition is justified because, rather than being driven by uterine volume changes, intrauterine pressure buildup is driven by increases in wall tension, as described by Laplace's law 38 : Since the length of every cell is given by Eq. (19) can be rewritten as We then have a nonlinear system of 3n unknowns and equations that needs to be solved at every time step. Solving this system, together with the cable equation, which includes the ionic current model of each cell in the myofibre, gives us the transmembrane voltage, the intracellular calcium concentrations, the length of each cell, and the force developed by each cell. Since all the model's variables are coupled, the entire system described by the equations in the "Methods" section needs to be solved at every time step to calculate the model outputs. For example, to calculate the length of every cell ( l (i) c ) from Eq. (20), the internal cells' lengths ( l (i) a , l (i) x , and l (i) s ) need to be calculated using Eqs. (15)(16)(17)(18)(19), which in turn require calculating the myosin phosphorylation states (Eqs. 6 and 7), and the cable Eq. (4), which requires solving the cells' ionic currents (Eqs. 1-3). Due to the complexity of this system, it is solved using the GNU scientific library in C++ 39 .
The excitation and contraction waves are coupled. While the electrical activity produces mechanical deformations and tensile force, the mechanical deformations feed back into the electrical activity through the geometrical parameters of the cable equation. This feedback is obtained through the left hand side of the cable equation, shown in Eq. (4).
To solve the cable Eq. (4) at every time step, we discretize the myofibre and consider each cell as an infinitesimal element. Using discrete operators, the left hand side of Eq. (4) becomes: where the superscript (i) denotes the ith cell in the myofibre. This expression closes the feedback loop since the length of the ith cell in the myofibre is used to calculate the electrical propagation across that cell and depends on the mechanical state of that cell.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.