An electro-mechanical multiscale model of uterine pregnancy contraction

Computers in Biology and Medicine 77 (2016) 182–194

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

An electro-mechanical multiscale model of uterine pregnancy contraction Maxime Yochum n, Jérémy Laforêt, Catherine Marque Sorbonne University,Université de technologie de Compiègne, CNRS, UMR 7338 Biomechanics and Bioengineering, Centre de recherche Royallieu, CS 6031960203 Compiègne cedex, France

art ic l e i nf o

a b s t r a c t

Article history: Received 5 February 2016 Received in revised form 30 July 2016 Accepted 1 August 2016

Detecting preterm labor as early as possible is important because tocolytic drugs are much more likely to delay preterm delivery if administered early. Having good information on the real risk of premature labor also leads to fewer women who do not need aggressive treatment for premature labor threat. Currently, one of the most promising ways to diagnose preterm labor threat is the analysis of the electrohysterogram (EHG). Its characteristics have been related to preterm labor risk but they have not proven to be sufficiently accurate to use in clinical routine. One of the reasons for this is that the physiology of the pregnant uterus is insufficiently understood. Models already exist in literature that simulate either the electrical or the mechanical component of the uterine smooth muscle. Few include both components in a co-simulation of electrical and mechanical aspects. A model that can represent realistically both the electrical and the mechanical behavior of the uterine muscle could be useful for better understanding the EHG and therefore for preterm labor detection. Processing the EHG considers only the electrical component of the uterus but the electrical activity does not seem to explain by itself the synchronization of the uterine muscle that occurs during labor and not at other times. Recent studies have demonstrated that the mechanical behavior of the uterine muscle seems to play an important role in uterus synchronization during labor. The aim of the proposed study is to link three different models of the uterine smooth muscle behavior by using co-simulation. The models go from the electrical activity generated at the cellular level to the mechanical force generated by the muscle and from there to the deformation of the tissue. The results show the feasibility of combining these three models to model a whole uterus contraction on 3D realistic uterus model. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Uterine model Multiscale Multi-physic Co-simulation

1. Introduction Premature birth is 9.6% of all births worldwide making preterm birth one of the world's largest public health problems [1]. The highest rates of preterm birth are in Africa and North America (11.9% and 10.6% of all births, respectively) [1]. Premature birth leads to high morbidity and mortality of newborns and can also lead to long-term adverse consequences for baby's and parents' health [2]. Children born at term have lower rates of cerebral palsy, sensory deficits, learning disabilities and respiratory illnesses when compared with children who are born before term, which results in enormous physical, psychological and economic costs [3,4]. To prevent premature labor efficiently, it is important to detect preterm labor symptoms as soon as possible via some biochemical or biophysical indicator. The electrical activity of the uterus is one of the most promising tools for the detection of n

Corresponding author. E-mail address: [email protected] (M. Yochum).

http://dx.doi.org/10.1016/j.compbiomed.2016.08.001 0010-4825/& 2016 Elsevier Ltd. All rights reserved.

preterm labor signs [5–7]. The electrohysterogram (EHG), which is the uterine electrical activity recorded externally on a women's abdomen is representative of uterine contraction as it is directly related to the electrical activity of the muscle which is the trigger of the mechanical contraction [8,9]. The prediction of a preterm labor by processing the EHG is possible from the 28th weeks of gestation [10], time when preterm labor can be treated effectively if detected soon enough or the baby prepared for premature birth using steroids. The uterine electrical activity detected on the abdomen's surface is generated at the cellular level. It passes through different kinds of tissue including muscle, fat and skin before reaching the abdominal electrode and being recorded as EHG. A correlation clearly exists between the electrical and the mechanical activities of the uterine muscle [11]. The mechanical effect of the uterine contraction depends on the cells excitability, as well as on the synchronization of the uterine cell activities, which occurs for whole cells in a short time during labor. This uterine synchronization could be related to the appearance of communicating junctions (Gap Junction) and by electrical diffusion

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

through the extracellular medium, permitting thus the excitability to propagate over the uterus [12–14]. But, the electrical propagation over the tissue does not on its own explain this uterine global synchronization present during labor [15]. Recent studies suggest that this uterine synchronization is a key factor in developing efficiency of uterine contractions [16]. A new hypothesis has been proposed that explains this uterine synchronization by means of a hydrodynamic-stretch activation mechanism that involves the mechanical effects of contractions (Intra Uterine Pressure IUP increase, tissue stretching) as trigger of the electrical activity through a mechano/electric coupling [15]. Other links exist between the EHG and the mechanical effect of the contraction. These links are due to the tissue deformation resulting from the force generated by the smooth muscle cells (SMCs). The muscle contraction and the resulting tissue deformation modify the position of the cells. As the cells move due to these contractions the cells that are “seen” by the electrodes positioned on the skin change during contractions. This displacement of the electrodes with respect to the uterine muscle surface can cause artifacts on the EHG signal [17] and introduce nonlinearities in the signal. Another aspect that could be of interest is that the deformation of the tissue modifies the distance between adjacent cells, which may influence the electrical propagation by changing the tissue conductivity locally [18–20]. Furthermore, knowing in detail the forces generated by each active cell and their resulting tissue deformation will permit to compute the Intra Uterine Pressure which is another measure that physicians use during active labor to evaluate the uterus contractile activity. However, IUP cannot be used during pregnancy due to its invasiveness [21]. It has been proven that the electrical and the mechanical activities of smooth muscles cells (SMCs) are linked together [22]. It is widely admitted that the contraction of SMCs is initiated by the intracellular calcium concentration, named [Ca2 +]i , which links to the contractile elements of the SMC (actin and myosin proteins) [22]. It has also been proven that mechanical stretching of the SMC can induce modifications in the electrical behavior of the cell, especially through stretch-activated channels [23,24]. The aim of this paper is to present a novel model of the uterine smooth muscle behavior by means of the co-simulation of different models representing its different electrical and mechanical aspects. The objective is then to demonstrate the feasibility to combine three different models, at the uterine cell and tissue levels, in order to make a co-simulation of the contractile process at whole uterus level. This co-simulation could then allow us to link the EHG signal and the mechanical effects of contraction. The two major novelties of this study are that we use an electrical model with biophysical meaningful parameters instead of a simple neuron model (Fitzhugh-Nagumo), and that we applied this model on the geometry of a pregnant human uterus obtained from experimental data. The use of this geometry gives a more realistic simulation of the model. The first model is an electrical one, dedicated to the cell/tissue level. It represents the ionic currents responsible for the uterine SMC activity, and its propagation over the tissue through electrical diffusion (simulating gap junctions). The second one is a force generation model at the cell level. It represents the mechanical force generated by a SMC due to the sliding of actin and myosin proteins. This force is related to the calcium concentration resulting from the calcium behavior modeled by the electrical model. The third one is a mechanical deformation model, which represents the modification of the uterine muscle shape due to the SMC forces generated by the second model. By combining these three models, the uterine behavior is modeled here in multi-scale (from the SMC to the whole organ) and with a multi-physics (chemical, electrical and mechanical) approach.

183

Simulation results are given for the behavior of: one cell, a 2D matrix of cells and a realistic 3D mesh of the uterine muscle, obtained from MRI images. Results show the feasibility of co-simulating uterine smooth muscle and its electrical, force and deformation aspects.

2. Materials and methods The challenge in this work is to use several existing models and to couple them in order to model the uterine muscle contraction in its different physiological and physical aspects: chemical, electrical and mechanical. In this section, we present the three models used for this study: the electrical cell model, the force cell model and the whole uterus deformation model. These models are combined to form a multi-scale and multi-physic co-simulation of the uterine smooth muscle. 2.1. Electrical model The electrical model used in this study is a reduction of the Rihana's model, a cellular model developed in our lab. All details can be found in [25,26]. The reduction of Rihana's original model (10 equations/cell) was a necessary step toward the modeling at the organ scale (whole uterine muscle) as all the equations in the cellular model for each cell have to be computed, at each sample time. Calculating this high number of equations is problematic when simulating several thousand of cells even for recent computer. The electrical model reduction, called Red3, uses only three variables [27]. This reduction is based on the change of the original ICa calcium current, which was the most time consuming in the original model, by the one proposed by Parthimos et al. [28]. The three variables of the simpler electrical model are described by the following equations:

⎞ dVm 1 ⎛ = ⎜ Istim − ICa − IK − IKCa − IL⎟, dt Cm ⎝ ⎠

(1)

h K − nK dnK = ∞ , τnK dt

(2)

d⎡⎣ Ca2 +⎤⎦

i

dt

⎛ ⎞ ⎜ ⎟ = fc ⎜ −αICa − kCa⎡⎣ Ca2 +⎤⎦ ⎟, i⎟ ⎜ ⎝ ⎠

(3)

where Vm is the transmembrane potential (initial condition is
50 mV), nK is the potassium activation variable (initial condition is 0.079257), KCa is the Calcium extraction factor and [Ca2 +]i the intracellular calcium concentration (initial condition is 0.1 mM). The ionic currents are ICa for the voltage dependent calcium channel current, IK for the voltage dependent potassium channel current, IKCa for the calcium dependent potassium channel current and Ileak for the leakage current. Equations for each current are then:

1

ICa = Jback + GCa(Vm − ECa) 1+

IK = GknK (Vm − EK ),

⎛ VCa − Vm ⎞ ⎜ ⎟ e⎝ R Ca ⎠

, (4)

(5)

184

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

Table 1 Parameters of the electrical model.

Table 2 Parameters of the mechanical HaiMurphy model from Maggio et al. experiment [32].

Variable

Value

Description

Unit

Gk

0.064

Potassium channels conductance

mS/cm2

GkCa

0.08

Potassium/Calcium channels conductance

mS/cm

GL

0.0055

Leak channels conductance

mS/cm2

κ

0.01

Half-point potassium concentration

μmol

fc α

0.4 4e
5

calcium influx probability current conservation factor

mol cm2/μC

kCa

0.1

Ca extraction factor

EL EK R


20
83 8.314

Leak nerst potential Potassium nerst potential gas constant

T F ⎡ Ca2 +⎤ ⎣ ⎦e

295 96.487 3

Temperature Faraday constant Extracellular calcium concentration

JK−1mol−1 K kC mol μmol

Jback

0.023

Background calcium current

μA/cm2

GCa

0.022

VOCC conductance

VCa RCa


20.07 5.97

Half-point of the VOCC activation sigmoid Maximum slope of the VOCC activation

mS/cm2 mV mV

IKCa = GkCa

⎡ Ca2 +⎤2 ⎣ ⎦i ⎡ Ca2 +⎤2 + κ 2 ⎣ ⎦i

ms−1) mV mV

1+

⎛⎜ 4.2 − Vm ⎞⎟ e⎝ 21.1 ⎠

(6)

(8)

, (9)

⎛ −Vm ⎞ ⎜ ⎟

τnK = 23.75e⎝ 72.15 ⎠ .

(10)

The constant for the electrical model are given in Table 1 The diffusion process, which occurs when the transmembrane potential spreads toward neighboring cells, is described by an electrical coupling computation using the gap junction model proposed by Koenigsberger et al. [29]. The gap junctions are known to increase during parturition [30] and could then be an important component in the labor triggering process. A coupling

M

k1

k7

AM

Mp

k2 k4 k5 k6

⎡ Ca 2 + ⎤ ⎣ 1/2MLCK ⎦ k2 k3 k4 k5 k6 k7 nM K

464.0758 nM 0.1399 s
1 14.4496 s
1 3.6124 s
1 k2 k1 0.1340 s
1 4.7135 5.0859

(11)

where Cm is the membrane capacitance, Ix represents the different ionic currents and D is the diffusion tensor. The latter is dependent on the distance dk that exists between neighboring cells. Therefore:

(7)

1

hK ∞ =

Value

∑ Ix dVm = − ▿·D▿(Vm), dt Cm

Dk =

⎛ ⎡ 2 +⎤ ⎞ RT ⎜ ⎣ Ca ⎦e ⎟ ln⎜ , 2F ⎜ ⎡⎣ Ca2 +⎤⎦ ⎟⎟ ⎝ i⎠

Parameter

term is thus added to Eq. (1):

(Vm − EK ),

IL = GL(Vm − EL ),

ECa =

2

k3

AMp

Fig. 1. Diagram of the cross-bridge cycle of Hai-Murphy model. M means unattached and unphosphorylated, Mp means unattached and phosphorylated, AMp means attached and phosphorylated, AM means attached and unphosphorylated. k1–k7 are the transition rates from one state to another.

1 , RCmdk2

(12)

where R is the membrane resistance and k a couple of neighboring cells [31]. When modeling a matrix of cell on a Cartesian grid, the distances dk are considered as being constant and equal. Then the values of D are the same for all the cells in the modeled tissue. When a realistic 3D uterine mesh is modeled, these distances are not identical anymore. On a mesh, the diffusion is computed for each edge of the mesh. The value Dk needs thus to be computed for every pair of adjacent cells for the whole mesh. Obviously, when a deformation is applied to the tissue, regardless if the tissue is a 2D Cartesian matrix or a more complex mesh, the distances used for the diffusion tensor has to be computed at each integration time of the deformation model. Dk increases when the distance between cells decreases and vice versa. Therefore, a feedback exists here from the tissue deformation toward the electrical tissue model. Concerning the stimulation Istim, a constant stimulation current (Istim ¼ 0.1175) was applied to cells which are considered activated during a certain period of time. In this way, the cell or group of cells will generate action potentials only during that period. 2.2. Force model We used in this study the modified Hai-Murphy model from Maggio et al. [32] who modified the original model [33] to better correspond to uterine SMCs. The Hai-Murphy model has already been used on several types of smooth muscle [34] for example smooth muscle in the airway [35]. The mechanical contraction is based on four-state kinetic cross-bridges where the fraction of cross-bridges in each state can be computed. Fig. 1 shows a diagram of the four state cross-bridges denoted M for unattached and unphosphorylated, Mp for unattached and phosphorylated, AMp for attached and phosphorylated, and AM for attached and unphosphorylated. k1–k7 are the transition rates from one state to another, where k1 and k6 are the myosin head phosphorylation rates, k2 and k5 are the dephosphorylation rates, k3 is the myosin head attachment rate, and k4 and k7 are the myosin head detachment rates. Table 2 gives the k1–k7 values used in this study, which correspond to the uterine SMC developed in reference [32]. From this reference k3, k4 and k7 are computed with functions determined relatively to the sliding of the actin and myosin

0 -20 -40 -60

F (mN) Fig. 2. Schematic implementation of the deformation model. Top right: Cartesian grid where muscular cells are organized on a 2D square grid. A Kelvin Voigt model is inserted between two neighboring cells (horizontally and vertically). Bottom right: mesh where a Kelvin Voigt model is inserted at each edge of the mesh. Bottom left: schematic representation of the Kelvin Voigt model.

⎛ dM ⎞ ⎛ ⎞⎛ ⎞ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ dt ⎟ ⎜ ⎟⎜ ⎟ ⎜ dMp ⎟ ⎜ −k1 ⎟⎜ M ⎟ k2 0 k7 ⎜ ⎟ ⎜ ⎟⎜ Mp ⎟ k4 0 ⎜ dt ⎟ = ⎜ k1 −k2 − k 3 ⎟⎜ ⎟, ⎜ dAMp ⎟ ⎜ 0 ⎟⎜ AMp⎟ −k 4 − k5 k3 k6 ⎜ ⎟ ⎜ −k6 − k7⎟⎟⎜⎜ AM ⎟⎟ 0 k5 ⎜ dt ⎟ ⎜ 0 ⎜ dAM ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎝ dt ⎠ ⎝ ⎠⎝ ⎠

where the initial conditions of the four states are M(0) ¼1 and Mp (0)¼ AMp(0)¼ AM(0) ¼0 and the sum of the four states is always one. The link between the two cellular models (electrical and force models) is mediated through the calcium concentration given by the electrical model resulting in a force acting on the mechanical model. This link has been introduced by Bursztyn et al. [36] on the k1 parameter (therefore also k6) which corresponds to the phosphorylated rate of cross-bridges, and permits the passage from the M to Mp states and from AM to AMp. k1 is given by:

k1(t ) =

⎡ Ca2 +⎤nM ⎣ ⎦i ⎡ Ca2 +⎤nM + ⎡ Ca 2 + ⎤nM ⎣ ⎦i ⎣ 1/2MLCK ⎦

0

5

10

15

5

10

15

20 c.

0.02 0

20 b.

0.04

0

5

10

Time (s)

15

20

the equation of the Kelvin Voigt model can be given as:

1 σ1, E

dε 1 = σ2, dt η

σ = σ1 + σ2,

(16)

where s1 is the stress of the spring, s2 is the dash-pot stress and s is the total model stress, E is the stiffness coefficient and η the viscosity of the material. The stretch between two cells is then computed by ε Eq. (17).

σ = Eε + η

(13)

0

Fig. 3. Example of the co-simulation of both electrical and force models when one AP is induced. (a) transmembrane potential Vm. (b) calcium concentration obtained from the electrical model. (c) force generated from the force model, related to the calcium concentration provided by the electrical model.

ε= proteins. However, in this study, the attachment and detachment rates k3, k4 and k7 are considered as constant and their values correspond to the average rates of attachment and detachment determined in [32]. In addition, k1¼k6 and k2¼ k5 and the matrix form equation of cross-bridges distribution is given by:

0

185

a.

-3 5 × 10

[Ca2+]i (mM)

Vm (mV)

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

dε dt

(17)

This model permits the connection between neighboring cells. The forces generated from the force model (presented Section 2.2), are connected to the Kelvin Voigt model through external forces (Fext) computed for each cell, according to the force generated by each given cell and the forces generated by its neighbors. The external force of a given cell is then obtained by a tensor computed as follow

Fext =

N

⎛

⎞⎛

⎞

i=1

⎝

⎠⎝

⎠

∑ ⎜⎜ Ci − C0⎟⎟·⎜⎜ F0 + Fi⎟⎟,

(18)

where, N is the number of neighbors for the concerned cell, C is the coordinate vector which can be a 2 values long vector in the 2D case, or a 3 values long vector in 3D case, and F is the force

, (14)

2+ where [Ca1/2 MLCK ] is the concentration of intracellular calcium that corresponds to the half activation of the myosin light-chain kinase, and nM is the Hill coefficient of this activation. The force generated by the SMC is the sum of both attached states (AM and AMp) with a proportional coefficient K.

F (t ) = K (AMp(t ) + AM (t ))

(15)

2.3. Deformation model In a simple first approach, we used a deformation model based on the Kelvin Voigt model, which is composed of a spring and a dash-pot. This model is quite often used to model viscoelastic biological materials [37]. A diagram of the Kelvin Voigt model is depicted Fig. 2 (bottom left). From the two elements of this model,

Fig. 4. Example of the co-simulation of both electrical and force models where a burst of APs is induced during 15 s (a) transmembrane potential Vm and (b) calcium concentration obtained from the electrical model. (c) Force generated from the mechanical model, related to the calcium concentration provided by the electrical model.

186

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

generated from the force model. Fext forces (in mN) are then applied to each corresponding vertex in order to induce contractions on the tissue. The Kelvin Voigt model takes as input the Fext force applied to each of the vertex in the three dimensions, which are applied to the stress s, and the current strain at a certain time (due to the previous deformation) to compute the new strain ε, which is used to give the new vertex position in three dimensions. The Kelvin Voigt tends then to bring back the tissue in its initial position when no more force is generated in the tissue. Fig. 2, top right, presents a diagram of the deformation model implemented for a Cartesian grid, where muscular cells are organized in a 2D square grid. A Kelvin Voigt model is thus inserted between two neighboring cells (horizontally and vertically). Fig. 2, bottom right, presents the implementation of the deformation model for a mesh. A Kelvin Voigt model is inserted at each edge of the mesh. The purpose of this paper is to show the feasibility of connecting an electrical, a force and a simple deformation model. Therefore the force generated by each cell, given by the force model, and calculated from the intracellular Calcium concentration (provided by the electrical model), is an input into each node in the Kelvin-Voigt deformation model. 2.3.1. Implementation information The model was implemented with Python 2.7 and use classical python scientific libraries provided by Anaconda distribution. The differential equations were solved using the Euler method with a time step of 0.5 ms. To reduce memory usage, the result of the model is not saved at the same frequency as the Euler resolution method. They are saved with a lower frequency, which can be selected by the user. For the result presented here, the frequency for saving data is 10 Hz except for the propagation velocity result which was saved at 2 kHz in order to keep a good enough time resolution. The source code of the Cartesian version of this model can be found at https://gitlab.utc.fr/jlaforet/uemg_squared. This version includes the co-simulation of models on a Cartesian grid in 1D, 2D and 3D. The co-simulation is set by means of Json files and is saved to hdf5 files. A graphical interface permits an easy visualization in 3D of the results. The model can be installed as a pure python package and the simulation can be run from the interpreter (requires some knowledge of python). We also provides linux/macOSX/Windows executables for users who just wish tu run the simulation without programming knowledge.

3. Results The previous sections introduce the 3 models used to represent the different components of the uterine smooth muscle. This section will present several results obtained from the co-simulation of these models. We first simulated one cell, to understand both electrical and mechanical phenomena at the cell level, then a matrix of cells to simulate the tissue level. Finally, we computed the 3 models on a realistic 3D mesh of the uterine muscle, in order to see the integration of electrical and mechanical effects of the contraction at the whole organ level. 3.1. Cellular level: one cell Fig. 3 presents an example of the co-simulation of both electrical and force models on a single uterine SMC. A single action potential AP (Fig. 3a) is induced in the electrical model with a stimulation current Istim ¼0.1175 during the first 0.6 s added to Vm in Eq. (1). We can see, by means of the electrical model, that the AP triggers an increase in the calcium concentration (Fig. 3b). The force then starts to increase and remains active even when the calcium concentration decreases and slowly returns to its initial value, giving rise to a force twitch (Fig. 3c). This is related to the dynamics of the mechanical, which is much slower than the electrical one. This phenomenon is particularly present in the uterine muscle because there is a very strong attachment between the actin and myosin proteins [38]. We note that the generation of force is maintained even after the full repolarization of the cell and in the total absence of an action potential. Fig. 4 presents an example of both electrical and force model co-simulation when applying, to a unique cell, a continuous stimulation that induces the generation of a burst of APs. Istim ¼0.1175 during the first 15 s Fig. 4a presents the transmembrane potential Vm with 20 APs generated over 15 s. Each AP triggers an increase in the calcium concentration (Fig. 4b). In the mechanical model, the presence of the successive APs creates an increase in the force generated due to the temporal summing of the individual twitches generated from each AP. Indeed, Fig. 4c shows the force generated becoming larger and larger after each AP. This phenomenon is due to the dynamics of the calcium concentration, which is present in the cell at a high value during a much longer time than the electrical depolarization, making the force generation lasts much longer. With continuous AP firing the force model reacts as if the cell is in tetanic contraction. This

Fig. 5. Co-simulation example with a two-dimension tissue (200  200 SMCs). Top: transmembrane potential Vm temporal signals of one cell (position: [50:50]). Bottom: six snapshots of Vm in 2D corresponding to time ¼ [0.05 s, 0.35 s, 0.65 s, 0.95 s, 1.25 s, 1.55 s].

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

187

Fig. 6. Co-simulation example with a two-dimension tissue (200  200 SMCs). Top: temporal force signals of one cell (position: [50:50]). Bottom: six snapshots of Force in 2D corresponding to time ¼ [0 s, 3 s, 6 s, 9 s, 12 s, 15 s].

cumulative effect on the force generation associated with a burst of AP is coherent with the bibliographic data. Indeed, an increase in the frequency of the APs, as well as of the EHG frequency content, has always been noticed when going from pregnancy to labor contractions and has been thought to reflect the increased strength of the contractions [8]. 3.2. Tissue level: a matrix of cells To simulate tissue behavior, we created a two dimensional tissue as a matrix of cells. Each point of the 2D matrix represents a uterine SMC, which is associated with an electrical and force model. Each cell is connected to its direct neighbors (top, bottom, left and right cells) with the simulation of gap junctions by the diffusion process (Section 2.1, Eq. (11)). First, the co-simulation is done without any deformation of the tissue, in order to simulate only the coupling between electrical activity of the cells and the resultant force. Then, the deformation model is added in order to evidence that the force, generated by the model, affects the position of cells. Even if the existence and localization of pacemaker cells is up for debate [39], we consider in this study as pacemaker area, one or multiple cell(s) that are electrically controlled by a stimulation current Istim in order to initiate action potentials in a chosen part of the tissue. 3.2.1. Without deformation Figs. 5 and 6 present an example of the co-simulation of both electrical (Fig. 5) and force (Fig. 6) models for a two-dimension tissue (200  200 SMCs). We control the stimulation of a pacemaker cell placed arbitrary at the center of the tissue (cell coordinates [100:100], Istim ¼0.1175 during the first 6 s). The temporal signals Vm (from the electrical model) and force (from the mechanical model) are plotted respectively in Fig. 5-top and Fig. 6top for the cell located at position [50:50]. We can notice the difference of dynamics between both models and the cumulative effect of the force generation. The bottom of Fig. 5 shows snapshots of the transmembrane voltage Vm obtained from the electrical model, corresponding to simulation times [0.05 s, 0.35 s, 0.65 s, 0.95 s, 1.25 s, 1.55 s], Note that the APs spread from the middle of the tissue toward the edges, due to the electrical diffusion between cells. The bottom of Fig. 6 shows snapshots of the force generated by the mechanical model. We can see that the force is rather homogeneous on the tissue (SMCs contracted on the whole tissue), while for Vm, the

potential value is very different between cells and at different times, as the AP spread out over the surface. This is caused by the slower dynamics of the mechanical model. This behavior is similar to the one of real uterine contractions, where a burst of EHG (containing lots of APs) is associated with a much slower single increase in intrauterine pressure [39]. A video of Figs. 5 and 6 is available online at: http://www.utc.fr/jlaforet/Suppl/2delecmeca. gif. 3.2.2. With deformation Fig. 7 presents an example of the co-simulation of electrical, force and deformation models for a two-dimension tissue (50  50 SMCs). We control the stimulation of a pacemaker cell located arbitrary at the center of the tissue (cell coordinate [25:25], Istim ¼0.1175 during all the simulation.). Fig. 7 presents 4 snapshots at different times of the co-simulation (t¼[0.2 s, 0.4 s, 0.6 s, 0.8 s]). To the right of each sub-figure the image at the top represents Vm (normalized), where we can see the APs spreading all over the tissue; the image at the bottom represents the force generated, where we can notice the presence of a force increment after the passing of each new AP. Fig. 7 (to the left of each subfigure) presents also the position of the 625 cells with their [x:y] coordinates (blue dots) corresponding to the lower left quarter of the simulation. Indeed, as the field is symmetric with respect to the central point of the tissue for this particular case, we represent only one quarter of the coordinates in order to be able to show more detail. The cell coordinates are computed by the deformation model (Kelvin Voigt model). In addition, we represent also quiver plots, depicted by the colored arrows, which represent the displacement of each cells. The base of each arrow represents the original coordinate of the cell and the arrow head represents the coordinate of each point at the simulation time. Due to the gradient of force (consequence of the propagation of the APs, which is a bit larger at the center of the tissue than at the edges), the distances between cells are smaller at the center of the tissue than at the edges. We can thus notice a contraction of the tissue. On this co-simulation, we keep the surface of the tissue constant by fixing the edges. A video of Fig. 7 is available online at: http://www.utc. fr/  jlaforet/Suppl/2delecmecadeformation2.gif. 3.3. Organ level: surfacic mesh of the uterus muscle We present in this section the co-simulation of the electrical and the 2 mechanical models on a realistic uterine mesh. The mesh was obtained from the FEMONUM project [40] (http://fe

188

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

Fig. 7. Co-simulation example with a two-dimension tissue (50  50 SMCs). The 4 snapshots are taken at different times of the co-simulation (a) t ¼0.2 s, (b) t ¼ 0.4 s (c) t¼ 0.6 s and (d) t¼ 0.8 s). To the right of each sub-figure Vm (normalized), where we can see the APs spreading all over the tissue (top), and the force generated. we can notice the presence of a force increment after each new generated AP (bottom). To the left of each sub-figure, the position of the 625 cells with their [x:y] coordinates (blue dots) corresponding to the lower left quarter of the simulation. Indeed, as it exists symmetry with respect to the central point of the tissue in this particular case, we represent only one quarter of the coordinates to not overload the figures. In addition, we represent also the quiver plots, depicted by the colored arrows, which represent the displacement of each cell. The base of each arrow represents the original coordinate of the cell and the arrow head represents the coordinate of each point for the current co-simulation time. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

monum.telecom-paristech.fr/) that offers the scientific community 3D fetal, uterine and abdominal meshes extracted from MRI imaging. Fig. 8 displays the original mesh including the fetus mesh (head downwards), seen from the mother's left side. The initial mesh contains 9288 vertices and 18,572 faces. However, this number of vertices seems too small to permit a proper simulation of the tissue. Indeed, when considering only 9277 vertices we get to a mean surface of 18.59 mm2 per vertex (cell). Previous work investigating the architecture of uterine smooth muscle [41] pointed out that the smooth muscle cells are organized in bundles, themselves organized in fasciculi. According to their measures, the fasciculi are 1–2 mm in diameter, which gives a section of 0.79–3.14 mm2, if we consider circular sections. Therefore, as we consider each vertex of the mesh as representing a fasciculus, we re-sampled the original mesh (by means of a Ball Pivoting algorithm) to obtain a vertex size approximately equal to the size of a fasciculus. The re-sampled mesh contains 99 084 vertices leading then to a mean surface of 1.74 mm2 per vertex, which corresponds to a fasciculus of 1.5 mm in diameter. In this

mesh, each vertex is associated with an electrical and a force models. Each edge is associated with an electrical diffusion process and to a Kelvin Voigt model, in order to be able to simulate the electrical propagation and the mesh deformation. The co-simulation is first done without any deformation of the tissue, in order to simulate only the electrical coupling between electrical activity of the cells and the resultant force. Then, the deformation model is added in order to evidence how the force generated by the model, affects the uterine shape. 3.3.1. Without deformation We first co-simulated only the electrical and force models, on the realistic uterine mesh. The electrical diffusion is performed through Eq. (11) which is applied to each edge of the mesh. A pacemaker area is defined on a mesh node located in the fundus of the uterus (Istim ¼ 0.1175 during the first 10 s). Fig. 9, displays snapshots of the simulation results at progressing simulation times. The transmembrane potential Vm is presented in the two upper lines (a–j) of Fig. 9. The action potentials spread along the

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

189

on the part where the force is the strongest. A video of Fig. 10 is available online at: http://www.utc.fr/  jlaforet/Suppl/model_ uterus3D_elec_meca_stretch.gif. In terms of computation time, one minute long co-simulation takes about 30 h on a laptop (Windows 8.1, Intel Core i3-2370M CPU @2.40 GHz and 6 Go DDR3 1333 MHz RAM). The model was implemented with Python 2.7 language, in a purely mono-thread fashion.

Fig. 8. original uterine and fetal meshes view from left mother's side.

mesh due to the electrical diffusion. The associated force generated from the mechanical model is presented in the two lower lines (k–t) of Fig. 9. The force spreads to the entire uterus and is present even when the cells are no longer electrically active (as depicted by Vm which is high only when an AP is present). A video of Fig. 9 is available online at: http://www.utc.fr/jlaforet/Suppl/ model_uterus3D_elec_meca.gif. 3.3.2. With deformation We then co-simulated the electrical and the two mechanical (force and deformation) models on the realistic uterine mesh. In these simulations, the values of the diffusion parameters Dk are computed at each integration step as the distances between cells change due to the deformation of the tissue. The electrical diffusion between the cells of a given part should be faster and stronger when this part of the tissue is contracted due to the lower intercell distances. A pacemaker area is defined on the mesh, located at the fundus of the uterus Istim ¼ 0.1175 during the first 10 s (original uterine and fetal meshes presented Fig. 8). In order to simulate a realistic contraction (volume nearly constant due to the presence of the fetus and of the amniotic fluid), the force is directed outside the uterus wherever the force is lower than the mean force of all cells, and to the inside wherever the force is greater than the mean force. Snapshots of the simulation results at progressing simulation times present the transmembrane potential Vm (Fig. 10 a–e), the associated force generated from the force model (Fig. 10f–j), the deformation (Fig. 10k–o), and the stretching of the uterus (Fig. 10p–t). The stretching is computed for the edges of the mesh as the difference between original edge lengths and edge lengths at a current time. A grid was superimposed to Fig. 10q–x, where each blue dot represents a cell, in order to better evidence the tissue deformation. We can notice here that the cells move closer

3.3.3. Propagation velocity The previous results evidenced a propagation of the electrical cell activity along the uterus. That propagation can be observed in real situations by means of multiple EHG measurements. Particularly when using a matrix of electrode, which allows the recording of several EHGs simultaneously. One of the analysis that can be done on EHGs during labor is the computation of the propagation velocity (PV) [42]. The difficulty of quantifying the electrical propagation in the uterus when using EHGs is to be sure that the PV is calculated in the perpendicular direction with respect to the propagation front. However, when using a model, as presented here, the direction of the propagation front can be evidenced. In order to compute the PV, two cells have been chosen on the front of the uterus (beside the woman's navel) as being visually perpendicular to the propagation front, separated by a distance equal to 47.55 mm (see the mesh representation, Fig. 11). Then 10 co-simulations have been performed with the models presented in this paper: 5 have been done including the deformation model, and 5 others without any deformation. Fig. 11-left presents the 5 different locations of a pacemaker area that have been chosen to induce a perpendicular front AP propagation at the electrode locations. We can see Fig. 11-right examples of results where panels a. (cell1) and b. (cell2) correspond to a simulation without deformation and panels c. (cell1) and d. (cell2) correspond to a stimulation with deformation. Then, the PV is computed by using the time difference of the local maxima of the transmembrane potential Vm extracted from both cells. We can identify about ten to eleven maxima per simulation. From these simulations, we obtain two sets of PV data, one set with deformation and another one without deformation. PV means computed without or with deformation are respectively 30.71 mm/s (0.71 standard deviation) without and 29.69 mm/s (0.12 standard deviation) with using the deformation model. These velocities are coherent with the real values computed from real EHG [42,43] that range between 6.6 and 138 mm/s with a mean of 21.5 mm/s [43] or 21.8 (76.8) mm/s [42]. Furthermore, we can notice a difference in PV depending on the deformation model is used or not. The propagation velocity decreases when using the deformation model (
3.33%) showing the influence of the increase in inter-cellular distances (at the position of the uterus which is observed to compute PV) linking to the diffusion tensor Eq. (11). We used ANOVA for a single factor to compare the two cases. The F obtained is 48.5, which corresponds to a p-Value equal to 2.1 × 10−10 . This indicates a statistically significant difference between the 2 the PV values computed with or without using the deformation model.

4. Discussion We presented in this paper a co-simulated model of the uterine muscle behavior developed with a multi-scale and multi-physics approach. The uterine muscle is modeled in its different aspects by using three models: one electrical (cell level), and two mechanical: force (cell level) and deformation (tissue/organ level). The final aim of this work is to better understand the synchronization

190

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

Fig. 9. Result of electrical and force models applied to a realistic 3D uterine mesh corresponding to the Fig. 8 with a view from left mother's side. Sub-figures (a) to (j) correspond to the electrical activity of the uterus and sub-figures (k)–(t) correspond to the mechanical activity of the uterus. A pacemaker area initiates electric action potentials at the fundus of the uterus.

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

191

Fig. 10. Result of electrical, force and deformation models applied to a realistic 3D uterine mesh corresponding to the Fig. 8 with a view from left mother's side. Sub-figures (a)–(e) correspond to the electrical activity of the uterus, sub-figures (f)–(j) correspond to the mechanical activity of the uterus, sub-figures (k)–(o) correspond to the mechanical deformation of the uterus (a grid was added for a better visualization), and sub-figures (p)–(t) correspond to the stretch of the tissue (blue is for contracted part and red is for relaxed part). A pacemaker area initiates electric action potentials at the fundus of the uterus. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

192

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

Fig. 11. Result of the propagation velocity computation. The position of the two considered cells is shown on the uterine muscle (left) as well as the location of the 5 pacemaker areas used for each simulation. The electrodes are parted away by 47.55 mm; a and b: Vm transmembrane potentials given for cell1 (respectively cell2) by the simulation without including the deformation model; c and d: Vm transmembrane potentials given for cell1 (respectively cell2) by the simulation including the deformation model. The red dots on these figures represent the local maxima used to compute the propagation velocity. Each red dot on a (and c) is related to the corresponding red dot in b (and d), permitting to compute the time difference between the 2 cells. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

process that occurs at the uterine level during labor. Previous studies allowed us to develop an electrical model of the uterine muscle cell behavior, which represents the ionic current dynamics leading to the cell contraction, as well as the electrical diffusion of the activity from cell to cell. However, when studying the uterine electrical activity (EHG) together with the global muscle contraction (IUP), Young stated that other aspects than the electrical diffusion alone must account for this synchronization process [15]. This paper demonstrates the feasibility of coupling existing models to co-simulate these electrical, force and deformation behaviors of the muscle, applying them to a 3D uterine muscle mesh representing a realistic uterus. The electrical model, adapted from the Hodgkin-Huxley model, describes the ionic currents involved in the uterine cell contraction, as well as the electrical diffusion of the activity to neighboring cells. It also provides the intracellular calcium concentration as an output parameter, which is then used as an input parameter to the force model. The force model gives, as output, the force generated by each active cell through the activation of actin and myosin protein sliding controlled by the intracellular calcium concentration. Then a simple deformation model (Kelvin Voigt) is used to simulate the deformation of the uterine muscle at the tissue and organ levels by moving the cells depending on the force they generate. The co-simulation of these 3 models has been performed at different scales including the cell, the tissue and the organ levels. At the cell level, we noticed a difference of dynamics between electrical and force model behaviors. The force dynamics is slower than the electrical one (see Figs. 3 and 4). This difference could be explained by the strong attachment of the actin and myosin proteins, which delays the sliding back of those proteins to their original position [38]. The other characteristic, evidenced by the co-simulation, is the cumulative effect on the force developed by a cell associated with a burst of action potentials. One action potential produces a force twitch. Repetitive action potentials, present in a burst, induce a strong increase in the maximal force, due

to the temporal summation of the overlapping force twitches. This result is very important and validates the principle of our uterine model. Indeed, many studies have evidenced a higher frequency content of the EHG during inefficient contractions (small increase in IUP), recorded during pregnancy, with respect to efficient contractions (strong increase in IUP), recorded during labor [8]. As the EHG represents the summation of all the cells electrically active during a contraction, this increase in its frequency content is representative of the increase in the firing rate of the uterine cells also seen at the cellular level [13]. The present simulation results permit thus to relate this increase in the frequency content of the EHG to the increase in IUP. Indeed, the IUP originates from the mechanical summation of the force generated by the active cells. If a cell contracts in a more efficient way, when its AP firing rate has to have increased, it thus appears logical that the IUP increase is stronger (the contraction more efficient) when the EHG signal shows higher frequency content. Then, a co-simulation of the electrical, force and deformation models was performed at a tissue level. In this simulation, cells are connected to each other with both electrical diffusion and a hyperelastic deformation model (Kelvin Voigt model). We noticed at this level that the tissue is homogeneously contracted even when the electrical transmembrane potential is heterogeneous (or inactive) in the tissue (see Figs. 5 and 6). This simulation shows that the mechanical properties of the tissue could explain the increase in IUP, because a piece of uterine tissue is able to be completely contracted, even if the tissue is no longer electrically active. This behavior cannot be obtained if we consider only the electrical properties of the tissue, which are localized in time and space. Indeed, when considering the propagation of action potentials, the strip of tissue electrically activated by an AP is very thin when compared to the strip of tissue mechanically active at the same time. Furthermore, the deformation model at the tissue level makes the cells generating force move closer to their neighbors. As the distances between neighboring cells are used in the electrical

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

diffusion (parameter Dk, Eq. (11)), the co-simulated model functions in closed loop. This could explain a part of the uterine synchronization, by facilitating the diffusion as the cells become closer. This also leads to a minor but consistent variation in propagation velocity of the APs. The last example presented here is the co-simulations of the electrical, force and deformation models at the uterine muscle scale (organ level). For computational reasons, it is obviously impossible to model the number of cells contained in a real uterine muscle (several billions of cells). Therefore, we used a 3D mesh extracted from pregnant woman's MRI. This mesh gives a realistic uterine shape where each vertex was considered as a muscular cell (associated to one electrical and one force model) and the edges of the mesh were associated to both electrical diffusion and Kelvin Voigt models. Because in this mesh the distances between cells vary during simulation, the spatial diffusion of the electrical model has been computed for each edge of the mesh and at each integration time. Due to the difference of dynamics between the electrical and mechanical models, the result shows a generation of force over the entire uterus. The maximal force is obtained from the mechanical model, only after the generation of several action potentials, due to the summation of their related force. We also computed, from the organ level co-simulation, the propagation velocities with and without including the deformation model. The result shows that the propagation velocity of the model is coherent with the propagation velocity that can be found in the literature. The difference of PV values, between using or not using the deformation model, shows the importance of considering the feedback loop from the deformation model toward the electrical diffusion coefficient. The contraction of the tissue induces changes in the distances between cells, and then, modifies the electrical diffusion. The model described in this paper presents steps beyond the state of the art when compared to the recently published models of uterine contractibility [44,45]. First, we use an electrical model representing the smooth muscle cell behavior, which has been specifically developed for uterine smooth muscle cell [25]. It is not a simpler neural oriented FitzHugh Nagumo action potential model as in [45–49]. The interest of using uterine smooth muscle cell model is that it directly gives the calcium concentration, which is further linked with the force generation model of Hai Murphy [32]. This permitted us to directly connect the electrical model to the force generation model without passing through a chemical model to get the calcium concentration as it is done in [45]. However, it is clear that we used a reduced physiological model for the electrical representation of uterine cells. This reduced model enables us to compute the model at organ level in a reasonable amount of time but then suffer of a not perfect fidelity to all ionic flow through the cell membrane as we can find in [50]. Secondly, the models presented in the literature use non-realistic uterus shape, quite often a simpler shape approximation as a semi symmetric ellipsoid. Our model is able to work with realistic 3D mesh obtained from MRI images of pregnant women as shown in the result section. In order to analyze preterm labor from the co-simulation of the model, a further step will be to convert the electrical activity given by this model at the organ level, into EHG signal simulated at the skin surface. This could be done by applying a spatial filtering in order to simulate the effect of the volume conductor (muscle, fat and skin layers present between the myometrium and the electrode. Thanks to this simulation of EHG signals, it could be possible to identify the value of the physiological parameters of the model from real EHGs recorded on a given patient. The value of these physiological parameters, if out of the normal range, should help the doctor in detecting earlier a preterm labor threat and in treating it properly (model aided diagnosis approach).

193

5. Limitation of the co-simulated model and perspectives In the results presented for the co-simulation on a realistic 3D mesh, we let the action potentials spread uniformly to the whole uterine muscle. This phenomenon is not consistent with the real uterine muscle structure that is neither homogeneous nor isotropic. Indeed, in a real uterine muscle, cells are organized as bundles surrounded by collagen fibers [41]. These bundle structures are associated with a better propagation within these bundles than between bundles [41], which clearly could disturb or even stop the diffusion for some parts of the tissue. We have to further include in our model some inhomogeneity of the diffusion parameters in order to represent a more realistic anatomy of the uterine muscle. For instance, we could add an anisotropic behavior on the uterus by changing the electrical diffusion according to the direction of a pair of cell. Then, we could differentiate the longitudinal and circular diffusion on the uterus by tuning this diffusion parameter. Another major influence on the propagation of the electrical activity is the number of gap junctions, which play an important role in the electrical diffusion for the uterine muscle. It has been shown that gap junctions evolve along pregnancy. Early in pregnancy, gap junctions are few but leading up to labor the number of gap junction increases abruptly [51]. This kind of physiological phenomenon can easily be included in the co-simulated models by for instance modifying the diffusion parameters in the electrical model. We have included in our co-simulations a coupling between the electrical model and the force model via the intracellular calcium concentration. Several studies in the recent literature have also evidenced the existence of stretch activated channels on SMC membranes. These channels induce ionic currents triggered by tissue stretching [23,24]. We therefore plan to include in our model a feedback loop, going from the deformation model toward the electrical model to represent the effect of these stretch activated channels. Considering the amniotic fluid and the fetus as incompressible, the contraction of the uterine muscle should be done at a constant volume. Therefore, a contraction of a part of a tissue should induce stretch to another part. This stretching will then induce new action potentials due to the activation of these stretch activated channels. This phenomenon could be of great importance for explaining uterine muscle synchronization. This is the central idea in a recent hypothesis proposed by Young [15] of how the synchronization of the uterus in labor happens. The deformation models used in this first part of our project for our co-simulation are simplistic. But we intended by this choice to test the feasibility of simulating electro-mechanical coupling to model the uterine behavior. In order to be more realistic, and to be able to simulate simultaneously the IUP and the abdominal EHG starting from the cellular activity, we plan in our future work to model the uterus, the intrauterine volume and the mother's abdominal wall using finite element modeling coupled to the electric-mechanical model presented in this paper. Acknowledgment This work was carried out and funded in the framework of the Labex MS2T. It was supported by the French Government, through the program “Investments for the future” managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02).

Appendix A. Supplementary data Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.compbiomed. 2016.08.001.

194

M. Yochum et al. / Computers in Biology and Medicine 77 (2016) 182–194

References [1] S. Beck, D. Wojdyla, L. Say, A.P. Betran, M. Merialdi, J.H. Requejo, C. Rubens, R. Menon, P.F. Van Look, The worldwide incidence of preterm birth: a systematic review of maternal mortality and morbidity, Bull. World Health Organ. 88 (1) (2010) 31–38. [2] C. Huddy, A. Johnson, P. Hope, Educational and behavioural problems in babies of 32–35 weeks gestation, Arch. Dis. Child.-Fetal Neonatal Ed. 85 (1) (2001) F23–F28. [3] S. Petrou, The economic consequences of preterm birth during the first 10 years of life, BJOG: Int. J. Obstet. Gynaecol. 112 (s1) (2005) 10–15. [4] S. Petrou, Z. Mehta, C. Hockley, P. Cook-Mozaffari, J. Henderson, M. Goldacre, The impact of preterm birth on hospital inpatient admissions and costs during the first 5 years of life, Pediatrics 112 (6) (2003) 1290–1297. [5] C. Marque, J.M. Duchene, S. Leclercq, G.S. Panczer, J. Chaumont, Uterine ehg processing for obstetrical monitorng, IEEE Trans. Biomed. Eng. 12 (1986) 1182–1187. [6] C. Marque, J. Duchene, Human abdominal ehg processing for uterine contraction monitoring, Biotechnol. (Read., Mass.) 11 (1989) 187. [7] G. Fele-Žorž, G. Kavšek, Ž. Novak-Antolič, F. Jager, A comparison of various linear and non-linear signal processing techniques to separate uterine emg records of term and pre-term delivery groups, Med. Biol. Eng. Comput. 46 (9) (2008) 911–922. [8] D. Devedeux, C. Marque, S. Mansour, G. Germain, J. Duchêne, Uterine electromyography:a critical review, Am. J. Obstet. Gynecol. 169 (6) (1993) 1636–1653. [9] C. Buhimschi, M.B. Boyle, G.R. Saade, R.E. Garfield, Uterine activity during pregnancy and labor assessed by simultaneous recordings from the myometrium and abdominal surface in the rat, Am. J. Obstet. Gynecol. 178 (4) (1998) 811–822. [10] H. Leman, C. Marque, J. Gondry, Use of the electrohysterogram signal for characterization of contractions during pregnancy, IEEE Trans. Biomed. Eng. 46 (10) (1999) 1222–1229. [11] M.J. Berridge, Smooth muscle cell calcium activation mechanisms, J. Physiol. 586 (21) (2008) 5047–5061. [12] Y. Abe, Effects of changing the ionic environment on passive and active membrane properties of pregnant rat uterus, J. Physiol. 214 (1) (1971) 173. [13] H. Kuriyama, H. Suzuki, Changes in electrical properties of rat myometrium during gestation and following hormonal treatments, J. Physiol. 260 (2) (1976) 315–333. [14] S. Sims, E.E. Daniel, R. Garfield, Improved electrical coupling in uterine smooth muscle is associated with increased numbers of gap junctions at parturition, J. Gen. Physiol. 80 (3) (1982) 353–375. [15] R.C. Young, Myocytes, myometrium, and uterine contractions, Ann. N.Y. Acad. Sci. 1101 (1) (2007) 72–84. [16] T.Y. Euliano, D. Marossero, M.T. Nguyen, N.R. Euliano, J. Principe, R.K. Edwards, Spatiotemporal electrohysterography patterns in normal and arrested labor, Am. J. Obstet. Gynecol. 200 (1) (2009), 54-e1. [17] C. Rabotti, H. de Lau, N. Haazen, G. Oei, M. Mischi, Ultrasound analysis of the uterine wall movement for improved electrohysterographic measurement and modeling, in: Engineering in Medicine and Biology Society (EMBC), 2013 35th Annual International Conference of the IEEE, IEEE, 2013, pp. 7436–7439. [18] C. Peracchia, A.F. Dulhunty, Low resistance junctions in crayfish. Structural changes with functional uncoupling, J. Cell Biol. 70 (2) (1976) 419–439. [19] R. Garfield, D. Merrett, A. Grover, Gap junction formation and regulation in myometrium, Am. J. Physiol.-Cell Physiol. 239 (5) (1980) C217–C228. [20] J.E. Saffitz, A.G. Kléber, Effects of mechanical forces and mediators of hypertrophy on remodeling of gap junctions in the heart, Circ. Res. 94 (5) (2004) 585–591. [21] L.F. Bastos, M.F. Lobo, W.L. van Meurs, D. Ayres-de Campos, An intrauterine pressure generator for educational simulation of labour and delivery, Med. Eng. Phys. 32 (7) (2010) 740–745. [22] R.C. Webb, Smooth muscle contraction and relaxation, Adv. Physiol. Educ. 27 (4) (2003) 201–206. [23] C. Guibert, T. Ducret, J.-P. Savineau, Voltage-independent calcium influx in smooth muscle, Prog. Biophys. Mol. Biol. 98 (1) (2008) 10–23. [24] F. Sachs, Stretch-activated ion channels:what are they? Physiology 25 (1) (2010) 50–56. [25] S. Rihana, J. Terrien, G. Germain, C. Marque, Mathematical modeling of electrical activity of uterine muscle cells, Med. Biol. Eng. Comput. 47 (6) (2009) 665–675. [26] S. Rihana, J. Santos, S. Mondie, C. Marque, Dynamical analysis of uterine cell electrical activity model, in: Engineering in Medicine and Biology Society, 2006. EMBS'06. 28th Annual International Conference of the IEEE, IEEE, 2006, pp. 4179–4182.

[27] J. Laforet, C. Rabotti, J. Terrien, M. Mischi, C. Marque, Toward a multiscale model of the uterine electrical activity, IEEE Trans. Biomed. Eng. 58 (12) (2011) 3487–3490. [28] D. Parthimos, D. Edwards, T. Griffith, Minimal model of arterial chaos generated by coupled intracellular and membrane Ca2 þ oscillators, Am. J. Physiol.Heart Circ. Physiol. 277 (3) (1999) H1119–H1144. [29] M. Koenigsberger, R. Sauser, M. Lamboley, J.-L. Bény, J.-J. Meister, Ca2 þ dynamics in a population of smooth muscle cells: modeling the recruitment and synchronization, Biophys. J. 87 (1) (2004) 92–104. [30] R. Garfield, S. Sims, E. Daniel, Gap junctions: their presence and necessity in myometrium during parturition, Science 198 (4320) (1977) 958–960. [31] S. Rihana, Modélisation de l′activitéélectrique utérine (Ph.D. thesis), Compiègne, 2008. [32] C.D. Maggio, S.R. Jennings, J.L. Robichaux, P.C. Stapor, J.M. Hyman, A modified hai murphy model of uterine smooth muscle contraction, Bull. Math. Biol. 74 (1) (2012) 143–158. [33] C.-M. Hai, R.A. Murphy, Cross-bridge phosphorylation and regulation of latch state in smooth muscle, Am. J. Physiol.-Cell Physiol. 254 (1) (1988) C99–C106. [34] J. StÅlhand, A. Klarbring, G.A. Holzapfel, Smooth muscle contraction: mechanochemical formulation for homogeneous finite strains, Prog. Biophys. Mol. Biol. 96 (1) (2008) 465–481. [35] B. Brook, O. Jensen, The role of contractile unit reorganization in force generation in airway smooth muscle, Math. Med. Biol. 31 (2) (2014) 99–124. [36] L. Bursztyn, O. Eytan, A.J. Jaffa, D. Elad, Mathematical model of excitationcontraction in a uterine smooth muscle cell, Am. J. Physiol.-Cell Physiol. 292 (5) (2007) C1816–C1829. [37] N. Özkaya, M. Nordin, D. Goldsheyder, D. Leger, Mechanical properties of biological tissues, in: Fundamentals of Biomechanics, Springer, 2012, pp. 221– 236. [38] C.-M. Hai, H.R. Kim, An expanded latch-bridge model of protein kinase c-mediated smooth muscle contraction, J. Appl. Physiol. 98 (4) (2005) 1356–1365. [39] R.E. Garfield, W.L. Maner, Physiology and electrical activity of uterine contractions, in: Seminars in Cell and Developmental Biology, Elsevier, vol. 18, 2007, pp. 289–295. [40] L. Bibin, J. Anquez, J.P. de la Plata Alcalde, T. Boubekeur, E.D. Angelini, I. Bloch, Whole-body pregnant woman modeling by digital geometry processing with detailed uterofetal unit based on medical images, IEEE Trans. Biomed. Eng. 57 (10) (2010) 2346–2358. [41] R.C. Young, R.O. Hession, Three-dimensional structure of the smooth muscle in the term-pregnant human uterus, Obstet. Gynecol. 93 (1) (1999) 94–99. [42] L. Lange, A. Vaeggemose, P. Kidmose, E. Mikkelsen, N. Uldbjerg, P. Johansen, Velocity and directionality of the electrohysterographic signal propagation, PLoS One 9 (1) (2014) e86775. [43] E. Mikkelsen, P. Johansen, A. Fuglsang-Frederiksen, N. Uldbjerg, Electrohysterography of labor contractions: propagation velocity and direction, Acta Obstet. Et. Gynecol. Scand. 92 (9) (2013) 1070–1078. [44] A.L. Cochran, Y. Gao, A model and simulation of uterine contractions, Math. Mech. Solids 20 (5) (2013) 540–564. [45] B. Sharifimajd, C.-J. Thore, J. StÅlhand, Simulating uterine contraction by using an electro-chemo-mechanical model, Biomech. Model. Mechanobiol. (2015) 1–14. [46] O. Aslanidi, J. Atia, A.P. Benson, H. van den Berg, A.M. Blanks, C. Choi, S. Gilbert, I. Goryanin, B. Hayes-Gill, A.V. Holden, et al., Towards a computational reconstruction of the electrodynamics of premature and full term human labour, Prog. Biophys. Mol. Biol. 107 (1) (2011) 183–192. [47] E. Pervolaraki, A.V. Holden, Spatiotemporal patterning of uterine excitation patterns in human labour, Biosystems 112 (2) (2013) 63–72. [48] R.E. Sheldon, C. Mashayamombe, S.-Q. Shi, R.E. Garfield, A. Shmygol, A. M. Blanks, H.A. van den Berg, Alterations in gap junction connexin43/connexin45 ratio mediate a transition from quiescence to excitation in a mathematical model of the myometrium, J. R. Soc. Interface 11 (101) (2014) 20140726. [49] R.E. Sheldon, M. Baghdadi, C. McCloskey, A.M. Blanks, A. Shmygol, H.A. van den Berg, Spatial heterogeneity enhances and modulates excitability in a mathematical model of the myometrium, J. R. Soc. Interface 10 (86) (2013) 20130458. [50] W.-C. Tong, C.Y. Choi, S. Karche, A.V. Holden, H. Zhang, M.J. Taggart, A computational model of the ionic currents, Ca2 þ dynamics and action potentials underlying contraction of isolated uterine smooth muscle, PLoS One 6 (4) (2011) e18685. [51] R. Garfield, T. Tabb, Role of cell-to-cell coupling in control of myometrial contractility and labor, Control of uterine contractility, CPI Llc., Florida, 1994, pp. 39–81.