Multiscale forward electromagnetic model of uterine contractions during pregnancy
 Patricio S La Rosa^{1},
 Hari Eswaran^{2},
 Hubert Preissl^{2, 3} and
 Arye Nehorai^{4}Email author
https://doi.org/10.1186/17566649124
© La Rosa et al.; licensee BioMed Central Ltd. 2012
Received: 5 July 2011
Accepted: 14 October 2012
Published: 5 November 2012
Abstract
Background
Analyzing and monitoring uterine contractions during pregnancy is relevant to the field of reproductive health assessment. Its clinical importance is grounded in the need to reliably predict the onset of labor at term and preterm. Preterm births can cause health problems or even be fatal for the fetus. Currently, there are no objective methods for consistently predicting the onset of labor based on sensing of the mechanical or electrophysiological aspects of uterine contractions. Therefore, modeling uterine contractions could help to better interpret such measurements and to develop more accurate methods for predicting labor. In this work, we develop a multiscale forward electromagnetic model of myometrial contractions during pregnancy. In particular, we introduce a model of myometrial current source densities and compute its magnetic field and action potential at the abdominal surface, using Maxwell’s equations and a fourcompartment volume conductor geometry. To model the current source density at the myometrium we use a bidomain approach. We consider a modified version of the FitzhughNagumo (FHN) equation for modeling ionic currents in each myocyte, assuming a plateautype transmembrane potential, and we incorporate the anisotropic nature of the uterus by designing conductivitytensor fields.
Results
We illustrate our modeling approach considering a spherical uterus and one pacemaker located in the fundus. We obtained a travelling transmembrane potential depolarizing from −56 mV to −16 mV and an average potential in the plateau area of −25 mV with a duration, before hyperpolarization, of 35 s, which is a good approximation with respect to the average recorded transmembrane potentials at term reported in the technical literature. Similarly, the percentage of myometrial cells contracting as a function of time had the same symmetric properties and duration as the intrauterine pressure waveforms of a pregnant human myometrium at term.
Conclusions
We introduced a multiscale modeling approach of uterine contractions which allows for incorporating electrophysiological and anatomical knowledge of the myometrium jointly. Our results are in good agreement with the values reported in the experimental technical literature, and these are potentially important as a tool for helping in the characterization of contractions and for predicting labor using magnetomyography (MMG) and electromyography (EMG).
Background
Modeling the myometrium contractility during pregnancy is of clinical importance, since it can aid in understanding the mechanism of labor, and, thus, it can help in monitoring the health of both the fetus and the mother. The occurrence of labor is in general accompanied by the appearance of periodic contractions which increase the intrauterine pressure to the point that cervix dilatation is manifested [1]. However, from clinical experiences, not all uterine contractions are efficient, i.e., lead to labor. Term labor is expected to occur after the 37th week of pregnancy, but in the last decade the incidence of preterm labor has increased significantly (12% of all births). Preterm birth can cause health problems or even be fatal for the fetus if it happens too early, and, at the same time, it imposes significant financial burdens on health care systems [2]. Therefore, it becomes critical to better understand the mechanism of bioreproduction which, as a consequence, would allow for the development of more effective forms of therapy that might help to predict labor and control the occurrence of labor. In the following we will discuss briefly the uterine microanatomy, previous contractions models, and our modeling approach.
Uterine microanatomy
The uterine microanatomy is consistent with action potential propagation [4]: (i) myocytes are densely packed within a bundle, (ii) bundles are contiguous within a fasciculus, and (iii) fasciculi are contiguous via communicating bridges formed with myocytes. In addition, the uterine changes during gestation are accompanied by the formation of gap junctions, which are one of the mechanisms for transmitting of contractile activity from cell to cell in a coordinated manner [1, 4]. The structure of the fasiculata within the uterus has not yet been well defined, but generally it makes the propagation of the action potential anisotropic [5, 6].
Uterine contraction models
Uterine contractions can be described by their mechanical and electrophysiological aspects. A mechanical contraction is manifested as a result of the excitation as well as the propagation of electrical activities in the uterine muscle, and appears in the form of an intrauterine pressure increase. Existing models approach the problem separately at the organ level [7–12] or at the cellular level [13–17].
At the organ level, the models presented in [7–11] focus on predicting the contractile forces that closely resemble clinical measurements of normal intrauterine pressure during contractions in labor, and, more recently, the action potential propagation. In [7], the authors assume that the uterus is a hollow ovoid formed by discrete contractile elements that propagate electrical impulses, generate tension, and have defined contracting and refractory periods. The envisioned mechanism for intercellular communication is based on action potential propagation, which is simulated by using a discrete state model for each cell. In [8], the authors revisit the model developed in [7] and perform multiple experiments with different uterine shapes, cell numbers, and initial distributions of active and resting cells. In [9], the author uses a discrete state model for combining action potential propagation and intercellular calcium wave propagation, two mechanisms of intercellular communication. However, in [7–9], mathematical and physical descriptions of the models are not provided. In [11], the trajectories of growth and myometrial tension at the onset of labor are modelled using a statistical modeling approach. The authors model the trajectories of uterine wall thickness, volume, and tension as functions of gestational age using a prolate ellipsoid method, intrauterine pressure results from the literature, and ultrasound measurements of the shape of the uterus collected on 320 subjects. Regarding the electrical activity of uterine contractions, a model of myometrial action potential propagation as measured by surface electrohysterography is proposed in [10]. The authors develop a myometrium skin conduction model consisting of four parrallel layers (myometrium,abdomen,fat, and skin), and model the action potential using a Gamma probability distribution. This model assumes that the electrical conducitivity of the myometrium is isotropic and that the abdominal curvature is negligible, thus making it suitable to model the electrical potential in a limited area of the abdominal surface. Recently, in [12], the authors model the electrical propagation of action potential using a 3dimensional myometrium for different initial conditions and stimuli. The authors use a reactiondifussion equation coupled with a FizhughNagumo model of ionic currents and consider a homogeneous isotropic myometrium.
At the cellular level, the models focus on predicting the changes of ionic concentrations in the intracellular and extracellular media during a contraction, and, as a consequence, on modeling the transmembrane potential evolution of a myocyte as a function of time. In [13, 14] a model is developed to simulate the complete process of a single myometrial smooth muscle contraction, which is initiated by depolarization. The model is based on the electrophysiological properties of a myocyte and on the cellular mechanisms that relate the rise in concentration levels of intracellular ion calcium Ca ^{2 +} to stress production. In [17], the authors model all known individual ionic currents of myometrial smooth muscle close to labor and combined them into a mathematical model of myometrial action potential generation. The model is shown to successfully mimic several recordings of spontaneous AP and force in uterine smooth muscle.
The notational convention adopted in this paper is as follows: italic font indicates a scalar quantity, as in a; lowercase boldface indicates a vector quantity, as in a, except for vector fields used in Maxwell’s equations such as electric field E, magnetic field B, and current density J; upper case italic indicates a matrix quantity, as in A. The matrix transpose is indicated by a superscript “^{ T }” as in A ^{ T }, and the identity matrix of size n×n is denoted I _{ n }. The set ${\mathbb{S}}^{n}$ denotes the vector space of symmetric n×nmatrices, and the spaces of nonnegative definite matrices and positive definite matrices are denoted by ${\mathbb{S}}_{+}^{n}$ and ${\mathbb{S}}_{++}^{n}$, respectively. The inner product and norm defined in the Euclidean space are denoted by 〈·,·〉 and ∥·∥, respectively.
Methods
In this section we discuss the electromagnetic source representation of uterine contractions. We will introduce a fourcompartment volume conductor model formed by an anisotropic bidomain myometrium and will present the models for the extrauterine electrical potential, magnetic field, and myometrial current source density. Finally, we will present a numerical example to illustrate our modeling approach.
Volume conductor model
Extrauterine magnetic field and electrical potential
Solving the forward electromagnetic problem of uterine contractions requires the computation of B(r,t) and ϕ(r,t) at $\partial \mathcal{A}$ using Eqs.(1) and (3), assuming that J _{s}(r,t)is known in $\mathcal{\mathcal{M}}$ and G(r) is known in all domains defined by the volume conductor geometry (see Figure 3).
The biological current sources J _{s}(r,t)in the myometrium are the transmembrane ionic fluxes due to concentration gradients, which flow across the surface membrane of the myocyte (smooth cells) from the extracellular medium into the intracellular medium and vice versa. The density of these ionic currents is also referred to as the impressed current density, since its origin is non electrical in nature, and it is the primary cause for the establishment of an electric field that induces secondary density currents in a conductive domain. We will model J _{s}(r,t) using a bidomain approach, which has proved to be a successful method to study electrophysiological activity in the myocardium [21, 22] and, more recently, in the uterus [27].
Myometrial currentsource density model
where ${G}_{\mathcal{\mathcal{M}}}^{\prime}=\left({G}_{\mathrm{e}}^{\prime}+{G}_{\mathrm{i}}^{\prime}\right)\in {\mathbb{S}}_{++}^{3}$ is the bulk myometrium conductivity tensor. Since spatial variations of v _{m}(r,t)depend on the local establishment of a transmembrane current density, j _{m}(r,t)≠0, we define the impressed currentdensity source as ${\mathit{J}}_{\mathrm{s}}(\mathit{r},t)={G}_{\mathrm{i}}^{\prime}\nabla {v}_{\mathrm{m}}(\mathit{r},t)$. Note that J _{s}(r,t)exists only when the spatial gradient exists, i.e., only in a region where the myometrium is undergoing depolarization (excitation) or repolarization.
This set of reactiondiffusion equations is also known as the bidomain equations [21, 22]. The diffusion part of the equations governs the spatial evolution of both the transmembrane and extracellular potentials, and the reaction part is given by the local ionic current cell dynamics. The solutions for v _{m}(r,t) and ϕ _{e}(r,t) depend on J _{ion}(r,t), J _{stim}(r,t), and the conductivity tensors, in addition to the boundary and the initial conditions. Since our goal is to model the propagation of electrical activity in the myometrium, we are interested in the class of travelingwave solutions of these equations whose waveform depends on J _{ion}(r,t) and whose initiation depends on J _{stim}(r,t). In what follows, we describe the models for both current densities, J _{ion}(r,t)and J _{stim}(r,t), and for the conductivity tensors, ${G}_{\mathrm{i}}^{\prime}$ and ${G}_{\mathrm{e}}^{\prime}$.
Ionic current model
where ε _{1}, ε _{2}, k, v _{1}, v _{2}, v _{3}, δ, γ, and βare model constants, and w (in V) is a state variable of the model. The parameter ε _{1} (in Ωm^{2}) controls the sharpness of the leading and trailing edges of the action potential waveform; the smaller ε _{1} is, the more vertical the edge is. Note that ε _{1}is a quantity of resistivity, therefore the smaller its value, the greater the permeability of the membrane to ionic flux. The parameter ε _{2} (in s^{−1}) controls the action potential duration; the smaller ε _{2} is, the longer it takes a cell to recover. The parameters v _{1}, v _{2}, v _{3} (in V), and k (in 1/V^{2}) control the range of v _{m}(r,t). Note that for a given set of k,v _{1}, v _{2}, and v _{3}, the ratios β/γ and δ/γ control the excitability threshold of the cell. The larger β/γis, the lower the excitability threshold that sets the cell dynamic to an oscillatory stable behavior between resting and exciting states is. Over a certain value, the cell dynamic becomes bistable, that is, if the cell starts from a resting potential, it changes to an excited state and remains there. On the other hand, a very negative β/γ value results in a permanent resting state. In the Results Section, presented further below, we select the model parameters using phasespace analysis and using the transmembrane potentials recorded from isolated human myometrial strips at term as a reference [13, 30]. This model does not consider explicitly the Ca ^{2 +} dynamics, and, moreover, it assumes that changes in the intra and extracellular ion concentrations are insignificant even after several depolarizations. However, its simplicity facilitates the modeling of the propagation of depolarization waves in large 2D and 3D domains.
Stimulus current model
where h _{ i }(r,t) is a spatiotemporal function with range in [0,1], ν _{ i } is the amplitude (in V), and N _{p}is the number of pacemaker areas. Intuitively, the former should modify the excitability of the cell at a certain instant of time based on the threshold value. In particular, our model assumes that the uterine myocyte can act as either a pacemaker or pacefollower, specifically, the spontaneous electrical behavior exhibited by the myometrium is an inherent property of the uterine myocyte (see [19] for more details.) Note that the size, duration, and intensity of the pacemaker area need to be chosen such that a stable traveling waveform solution to the bidomain equations on the myometrium is possible.
Conductivity tensor model
The structure of the fasiculata within the uterus has not been completely characterized as a whole. Despite of not having a sense of a global uterinefiber arcuitecture, it has been possible to characterize local structures in each of the three layers of the uterine wall [5, 6, 12]. In [5], the authors investigated the global fiber architecture of the nonpregnant uterus by diffusion tensor magnetic resonance imaging (DTMRI). From the exvivo analysis of five nonpregnant uteri, the authors identified an inner circular layer around the uterine cavity on slices orthogonal to the long axis of the organ. In the regions outside the inner circular layer, they could not identify a global structure but did find several locally aligned groups of fibers. At the level of the cervix, they found an outer circular layer and an inner region with mostly longitudinal components. In [12], from the analysis of the DTMRI of one postpartum uterus, the authors concluded that the myometrium is organized into bundles, with the fibre bundles forming interweaving sheets. In the following, we will introduce an approach for designing the conductivity tensors in the myometrium.
Hence, to construct the conductivity tensors as a function of r, it is enough to define the vector field a _{3}(r), the myometrial fiber orientation, in each location of the anisotropic domain, as well as the conductivity values σ _{il}, σ _{el},σ _{it},and σ _{et}, because of the assumption of cylindric symmetry.
Designing myometrial fiber orientations
where $R=\sqrt{{x}^{2}+{y}^{2}+{z}^{2}.}$ Note that, for α=0, the main axis of the fibers runs vertically from the fundus to the cervix.
Estimating myocyte fiber conductivities
where p is the fraction of the volume occupied by the myocytes and collagenous fibers in the tissue, and m is the socalled cementation factor, which depends on the shape and orientation of the myocyte in the tissue.
Monodomain approximation and boundary conditions
The above simplification is also known as the monodomain approximation of the bidomain equations, which, under suitable boundary conditions, allows the computation of v _{m} and, thus, J _{s}, independent from ϕ _{e}.
where ${\hat{\mathit{n}}}_{j}$ is the normal vector to the surface j in each case.
Summary of modeling assumptions
In this section we list all the assumptions that were considered while developing our model:

Volume conductor geometry: it is given by four spherical compartments. The uterine wall is a single layerwith uniform thickness. The abdominal, intrauterine, and fetal compartments have isotropic and homogeneous conductivity. The uterine compartment is anisotropic and homogeneous. All boundaries between compartments, except for the external boundary of the abdominal and fetal compartments, are electrically conductive. The fetal compartment boundary can be set to be electrically insulated or electrically conductive to take into account the presence of vernix caseosa.

Myometrium: It is a continuous medium formed by array of myocytes with same transmembrane potential shape across the whole domain. The extracellular and intracellular tissue conductivities are anisotropic, homogeneous, and proportional to each other(equal anisotropic ratio). The fiber orientation with respect to the main symmetry axis is kept fix across the whole domain.

Myocyte: It is a cell of cylindrical shape. Its electrical conductivity is assumed to have a cylindrical symmetry with higher conductivity along the main axis of the cylinder than along the axis defining the cross section. Transmembrane potential shape is plateautype, and any myocyte can play the roll of pacemaker or follower by means of a stimulus current.
Numerical computation of v _{m}(r,t), ϕ(r,t), and B(r,t)
The computations of v _{m}(r,t), ϕ(r,t), and B(r,t)are given by the following procedure:

Step 1: Solve for v _{m}(r,t)using Eqs. (37), (12), (13), and (14) subject to boundary conditions (41) and (44), and to initial conditions given by$\begin{array}{l}\phantom{\rule{1em}{0ex}}{v}_{\mathrm{m}}(r,0)={v}_{\text{mr}},\hfill \end{array}$(47)$\begin{array}{l}\phantom{\rule{0.75em}{0ex}}{w}_{\mathrm{m}}(r,0)=\left(\frac{\beta}{\gamma}{v}_{\text{mr}}+\frac{\delta}{\gamma}\right),\hfill \end{array}$(48)$\frac{\partial {v}_{\mathrm{m}}(r,0)}{\mathrm{\partial t}}=0,\text{and}$(49)$\frac{\partial {w}_{\mathrm{m}}(r,0)}{\mathrm{\partial t}}=0.$(50)

Step 2: Solve for ϕ _{e}(r,t)in $\mathcal{\mathcal{M}}$ and ϕ(r,t) in $\mathcal{A}$ and $\mathcal{U},$ using the solution of v _{m}(r,t), computed in Step 1, and the following expressions:$\nabla \xb7\frac{\left(\varsigma +1\right)}{\varsigma}{G}_{\mathrm{e}}^{\prime}\nabla {\varphi}_{\mathrm{e}}(r,t)=\nabla \xb7{G}_{\mathrm{e}}^{\prime}\nabla {v}_{\mathrm{m}}(r,t),\phantom{\rule{2.77695pt}{0ex}}\text{in}\phantom{\rule{2.77695pt}{0ex}}\mathcal{\mathcal{M}},$(51)$\begin{array}{l}\phantom{\rule{3.5em}{0ex}}\nabla \xb7{G}_{\mathcal{A}}\nabla \varphi (r,t)=0,\phantom{\rule{2.77695pt}{0ex}}\text{in}\phantom{\rule{2.77695pt}{0ex}}\mathcal{A},\hfill \end{array}$(52)$\begin{array}{l}\phantom{\rule{3.5em}{0ex}}\nabla \xb7{G}_{\mathcal{U}}^{\prime}\nabla \varphi (r,t)=0,\phantom{\rule{2.77695pt}{0ex}}\text{in}\phantom{\rule{2.77695pt}{0ex}}\mathcal{U},\hfill \end{array}$(53)
subject to the boundary conditions (39), (40), (42), (43), (45), and (46).

Step 3: Solve for B(r,t)using Eq.(1), and computing the total current density J(r,t) in the whole uterine domain, using the solutions of v _{m}(r,t), ϕ _{e}(r,t), and ϕ(r,t), obtained in Steps 1 and 2, and considering the following boundary conditions:$\begin{array}{l}\phantom{\rule{1em}{0ex}}{\hat{\mathit{n}}}_{\mathcal{F}}\times \left({\mathit{B}}_{\mathcal{F}}(\mathit{r},t){\mathit{B}}_{\mathcal{U}}(\mathit{r},t)\right)=0,\text{in}\phantom{\rule{2.77695pt}{0ex}}\partial \mathcal{F},\hfill \end{array}$(54)${\hat{\mathit{n}}}_{\mathcal{U}}\times \left({\mathit{B}}_{\mathcal{\mathcal{M}}}(\mathit{r},t){\mathit{B}}_{\mathcal{U}}(\mathit{r},t)\right)=0,\text{in}\phantom{\rule{2.77695pt}{0ex}}\partial \mathcal{U},$(55)${\hat{\mathit{n}}}_{\mathcal{\mathcal{M}}}\times \left({\mathit{B}}_{\mathcal{A}}(\mathit{r},t){\mathit{B}}_{\mathcal{\mathcal{M}}}(\mathit{r},t)\right)=0,\text{in}\phantom{\rule{2.77695pt}{0ex}}\partial \mathcal{\mathcal{M}},$(56)$\begin{array}{l}\phantom{\rule{1em}{0ex}}{\hat{\mathit{n}}}_{\mathcal{A}}\times \left({\mathit{B}}_{\mathcal{A}}(\mathit{r},t){\mathit{B}}_{\mathcal{E}}(\mathit{r},t)\right)=0,\text{in}\phantom{\rule{2.77695pt}{0ex}}\partial \mathcal{A},\hfill \end{array}$(57)$\begin{array}{l}\phantom{\rule{7em}{0ex}}{\hat{\mathit{n}}}_{\mathcal{A}}\times {\mathit{B}}_{\mathcal{E}}(\mathit{r},t)=0,\text{in}\phantom{\rule{2.77695pt}{0ex}}\partial \mathcal{E},\hfill \end{array}$(58)
where $\mathcal{E}$ is an additional volume that must be incorporated in order to compute the magnetic field at the abdominal surface using FEM. The boundary conditions at the interfaces $\partial \mathcal{F},\partial \mathcal{U},\partial \mathcal{\mathcal{M}}$, and $\partial \mathcal{A}$ describe the continuity of the magnetic field between domains, and the boundary condition at $\partial \mathcal{E}$, the external boundary of $\mathcal{E}$, is of electric isolation.
To compute the solution in each of the above steps, we use the FEM solver COMSOL Multiphysics, running on a server with eight 64bit processors at 2.3 GHz, with 32 GB RAM.
Results
Conductivity values of the volume conductor geometry
Myocyte dimensions and Archie’s law parameters
In [19], the authors reported that, in the human uterus, there may be a preferential direction of propagation of contractions, and thus of transmembrane potential propagation, from the fundus toward the isthmus, which could aid in the expulsion of the fetus. Therefore, in order to study this assumption with our model, we consider j _{stim} with N _{ p }=1, ν _{1}=2 V, and h _{1}(r,t)={1, if 0≤t≤0.1, 0.15 ≤∥r∥ ≤0.16 and z≥0.15 ; 0,otherwise. The size and intensity of the pacemaker area are chosen in order to obtain a stable traveling waveform solution to the bidomain equations on the spherical myometrium.
To solve the set of equations that model the myometrial currentsource densities and their electromagnetic fields at the level of the abdomen, we use a threestep procedure along with Finite Element Methods (FEM) (See the Methods Section for more details). The FEM discretization of the whole volume conductor is done using tetrahedral elements. The elements used to discretize each domain consists of 75,964 elements divided as follows: 1,650 elements in the fetus, 7,605 elements in the intrauterine cavity, 8,299 elements in the myometrium, and 19,036elements in the abdominal cavity. To compute the magnetic field using FEM, an additional domain (domain $\mathcal{E}$ in (58)) concentric to the volume conductor geometry has to be added. Specifically, we add a concentric sphere with radius of 1m. The discretization of this domain consists of 39,374 elements. The discretized boundaries consists of 10,996triangular elements, and the number of degrees of freedom of the whole domain is 159,130. In all the Steps of the numerical computation procedure described above we use the generalized minimal residual method (GMRES) for solving the resulting system of linear equations after FEM discretization. The backward differentiation formula is used to discretize the time derivative in Step 1.
In Figure 9(b) we illustrate the percentage of the contracting myometrial volume as function of time, which in [7, 9] was used as a reference to compute the changes in the intrauterine pressure due to a contraction. Interestingly, we observe that the percentage of the contracting myometrial volume as a function of time has the symmetric properties and time duration of the intrauterine pressure waveforms of a pregnant human myometrium at term, as discussed in [9].
Discussion
Our numerical simulation results (see Figure 8) show that, because of the anisotropy in the conductivity, the direction of the current density J _{s}is rotated at a certain angle from the main direction of the propagating transmembrane potential v _{m}, and it is the transversal component of this current, parallel to the xy plane, which generates the magnetic field ${\mathit{B}}_{{\mathit{n}}_{\mathcal{A}}}.$ This observation is in agreement with the analysis presented in [37], and it is important to take it into account when interpreting the magnetic field measurements generated by uterine contractions in the presence of volume conductor geometry. Therefore, the spatial signature of ${\mathit{B}}_{{\mathit{n}}_{\mathcal{A}}}$ is highly dependent on the fiber orientation of the myometrium. Because of the proximity of the sensors and the myometrium, it is not strictly applicable to assume a moving dipole parallel to the direction of propagation of the transmembrane potential as the main model for the current source generating the measured magnetic field. This last interpretation might be suitable in the case where the transverse length of the transmembrane potential front is short in comparison to the area covered by the array of sensors. In contrast, if the transverse length is larger and thus not covered by the measured area, (for example, when several cells are recruited), then it is suitable to consider a moving line source (stretched ring) model instead.
Though our volume conductor geometry is an oversimplification of a real anatomical structure, certain aspects of it stand as a fair approximation of the volume conductor geometry when performing MMG measurements. For example, the SARA system, known as the superconducting quantum interference device array for reproductive assessment ^{a}(SQUID Array for Reproductive Assessment) is a passive, stationary, floormounted instrument at which the patient sits and leans her abdomen against a concave surface which contains an array of 151 sensors. The effect of leaning the abdomen on the concave surface works as a way to standardize the abdominal surface making the spherical model very suitable to represent the measuring surface. Also, a spherical uterus is good approximation at the early stage of pregnancy but not necessarily through all stages. In this sense, a more realistic populationaverage uterine model should be constructed from magnetic resonance images. Unfortunately, exposing of pregnant patients to MRI is not currently the norm, and thus having access to average uterine geometry for different stages of pregnancy is not yet possible.
The ionic current model based on extended FHN equations can reproduce a good approximation of the plateautype transmembrane potential recorded in human myocytes, however, it introduces hyperpolarization, which constrains the excitability of the cell, and thus consecutive contractions can only take place until v _{m} reaches resting potential. In our case, our model can simulate a minimum period of one contraction every 240 s, allowing us to model labor contractions whose period reduces to about one contraction every 300 s. In addition, note that a larger ε _{2} value can extend the transmembrane potential duration to values closer to the average duration reported in [13, 30]. However, a larger ε _{2}value also extends the duration of hyperpolarization and therefore the plateau of the curve describing the percentage of contracting myometrial volume. This last observation is the result of the combined effect of the transmembrane potential duration and the geometry of the uterus. Clinical observations on the shape of the uterine pressure wave have found correlation between the rising slope of the waveform and contraction efficiency [38], namely, the steeper the slope the more efficient the contraction is. Additional numerical examples, using one pacemaker at the same position, show that our model is consistent with the latter observation, since we found that changing the propagation speed of the activity modifies the shape of the percentage of the contracting myometrial volume as a function of time. In particular, a faster speed of propagation of the main transmembrane potential front increases the rising slope the percentage of the contracting myometrial volume, however, above a certain threshold it reduces the duration of the contraction. Placing two pacemakers at each extreme of the spherical domain, one at the fundus and the other one at the cervix, doubles the rising slope of the percentage of the contracting myometrial volume curve, and, depending on the duration of the transmembrane potential, it introduces a plateau in the curve, which implies that all volume is contracting. The symmetry of the contracting myometrial volume curve generated using our model is primarily due to the symmetry of the volume conductor used in our simulations and the duration of the transmembrane potential simulated.
Modeling the different stages of pregnancies can be accommodated using this model. First the volume conductor geometry should be scaled to the average volume size corresponding to stage of pregnancy of interest. Specifically, the fetal, the uterine, and the abdominal shapes should be modified accordingly using reported values in the literature (e.g., see [3] for uterine shapes). Second, the conductivity values of the volume conductor model should be set with respect to the stage of pregnancy. For example, the conductivity of the intrauterine cavity given by the amniotic fluid is known to change as the pregnancy develops (e.g., see [36]). Regarding the equal anisotropy ratio, this can be computed by using the reported values of the speed of propagation, the resting transmembrane potential, and the updated ionic current parameters such these fit the transmembrane potential for the specific stage. Third, the boundary conditions remain the same across all periods, noting that for boundary condition given by Eq. (46), λ should be set to 1 or 0 depending on the presence of vernix caseosa. Fourth, incorporation of a bursting type transmembrane potential can be done using the FHN model by adding a periodic function to Eq. (13), which control the excitability of the cell. However, more realistic and complex shapes of transmembrane potential can be designed by incorporating, for example, the ionic current model presented in [17].
Additional approaches to validate our mathematical model should consider comparisons against real magnetic and potential field measurements at the level of the abdomen, which leads automatically to the next natural question: how to solve the inverse problem of our model? Explicitly, given a set of measurements and a parametric model, can we infer the value of certain parameters of interest, such as the number of stimui and their positions, the conductivity tensor in the domain, initial and boundary conditions, etc? Solving the inverse problem stands as the connecting bridge between multiscale modeling of the pregnant uterus and clinical applications. In particular, as a first approach we are planning to use MMG measurements to estimate the current density in the myometrium and thus search for features that characterize contractions patterns which lead to pretermlabor.
Conclusion
We proposed a multiscaleforward electromagnetic model of uterine contractions during pregnancy. Our model incorporates knowledge of the electrophysiological aspects of uterine contractions during pregnancy at both the cellular and organ levels. We applied a bidomain approach for modeling the propagation of the myometrium transmembrane potential v _{m} on the uterus and used this to compute the action potential ϕand the magnetic field Bat the abdominal surface. We introduced a modified version of the FitzHughNagumo equation for modeling the ionic currents in each cell. Though our ionic current model does not consider Ca ^{2 +} dynamics explicitly, the simplicity of the modified equation allows for the propagating action potential to be modeled under well defined conditions as shown in this paper. We also proposed a general approach to design conductivity tensors in the myometrium and to estimate the conductivity tensor values in the extracellular and intracellular domains. We introduced a simplified geometry for the problem and proposed a discretized model solution based on a finite element method approach. Finally, we illustrated our modeling approach through a numerical example by modeling uterine contractions at term. Our model is potentially important as a tool for helping characterize contractions and for predicting labor using MMG and EMG.
As part of our future work, we will investigate pearshaped uterine domains as a way to approximate better the uterine geometry. We will also include more accurate ionic current models as in [13–15, 17] and will consider spatial variations of the fiber orientation.
Endnote
^{a}SARA was built in collaboration with VSMMedTech Ltd., Canada and is installed at the University of Arkansas for Medical Sciences (UAMS) Hospital.
Declarations
Acknowledgements
The study was funded in part by NIH Grant NIBIB/1R01 EB00726401A2, National Science Foundation, Grant No. CCF0963742, and by Grant DFG BI 19550.
Authors’ Affiliations
References
 Chard T, Grudzinskas JG: The Uterus. 1995, Cambridge: University PressGoogle Scholar
 Lamont R: Looking to the future. BJOG: Int J Obstet Gynaecol. 2003, 110 (20): 131135. 10.1046/j.14710528.2003.00065.x.View ArticleGoogle Scholar
 Degani S, Leibovitz Z, Shapiro I, Gonen R, Ohel G: Myometrial thickness in pregnancy: longitudinal sonographic study. J Ultrasound Med. 1998, 17 (10): 661665.PubMedGoogle Scholar
 Young R, Hession RO: Threedimensional structure of the smooth muscle in the termpregnant human uterus. Obstet Gynecol. 1999, 93: 9499. 10.1016/S00297844(98)003457.PubMedGoogle Scholar
 Weiss S, Jaermann T, Schmid P, Staempli P, Niederer P, Caduff R, Bajka M: Threedimensional fiber architecture of the nonpregnant human uterus determined Ex Vivo using magnetic resonance diffusion tensor imaging. Anat Rec Part A, Discov Mol, Cell, Evol Biol. 2006, 288: 8490.View ArticleGoogle Scholar
 Young RC: Myocytes, myometrium, and uterine contractions. Ann New York Acad Sci. 2007, 1101: 7284. 10.1196/annals.1389.038.View ArticleGoogle Scholar
 Andersen HF, Barclay ML: A computer model of uterine contractions based on discrete contractile elements. Obstet Gynecol. 1995, 86: 108111. 10.1016/00297844(95)001114.View ArticlePubMedGoogle Scholar
 Barclay M, Andersen H, Simon C: Emergent behaviors in a deterministic model of the human uterus. Reprod Sci. 2010, 17 (10): 948954. 10.1177/1933719110376544.View ArticlePubMedGoogle Scholar
 Young R: A computer model of uterine contractions based on action potential propagation and intercellular calcium waves. Obstet Gynecol. 1997, 89: 604608. 10.1016/S00297844(96)005029.View ArticlePubMedGoogle Scholar
 Rabotti C, Mischi M, Beulen L, Oei G, Bergmans J: Modeling and identification of the electrohysterographic volume conductor by highdensity electrodes. Biomed Eng, IEEE Transac. 2010, 57 (3): 519527.View ArticleGoogle Scholar
 Sokolowski P, Giles W, McGrath S, Smith D, Smith J, Smith R: Human uterine wall tension trajectories and the onset of parturition. PLoS ONE. 2010, 5 (6): e1103710.1371/journal.pone.0011037.View ArticlePubMedPubMed CentralGoogle Scholar
 Aslanidi O, Atia J, Benson A, van den Berg H, Blanks A, Choi C, Gilbert S, Goryanin I, HayesGill B, Holden A, Li P, Norman J, Shmygol A, Simpson N, Taggart M, Tong W, Zhang H: Towards a computational reconstruction of the electrodynamics of premature and full term human labour. Prog Biophys Mol Biol. 2011, 107: 183192. 10.1016/j.pbiomolbio.2011.07.004.View ArticlePubMedGoogle Scholar
 Bursztyn L, Eytan O, Jaffa AJ, Elad D: Modeling myometrial smooth muscle contraction. Ann New York Acad Sci. 2007, 1101: 110138. 10.1196/annals.1389.025.View ArticleGoogle Scholar
 Bursztyn L, Eytan O, Jaffa AJ, Elad D: Mathematical model of excitationcontraction in a uterine smooth muscle cell. Am J Physiol  Cell Physiol. 2007, 292 (5): C1816—C1829PubMedGoogle Scholar
 Rihana S, Marque C: Modelling the electrical activity of a uterine cell, a mathematical model approach. Proc. The 3rd European Medical and Biological Engineering Conference. 2005, PragueGoogle Scholar
 Rihana S, Terrien J, Germain G, Marque C: Electrophysiological model of the uterine electrical activity. Med Biol Eng Comput. 2009, 47 (6): 665675. 10.1007/s1151700904334.View ArticlePubMedGoogle Scholar
 Tong W, Choi C, Karche S, Holden A, Zhang H, Taggart M: A computational model of the ionic currents, Ca2+ dynamics and action potentials underlying contraction of isolated uterine smooth muscle. PloS One. 2011, 6 (4): e1868510.1371/journal.pone.0018685.View ArticlePubMedPubMed CentralGoogle Scholar
 Eswaran H, Preissl H, Wilson JD, Murphy P, Robinson S, Lowery C: First magnetomyographic recordings of uterine activity with spatialtemporal information with 151channel sensor array. Am J Obstet Gynecol. 2002, 187: 145151. 10.1067/mob.2002.123031.View ArticlePubMedGoogle Scholar
 Garfield RE, Maner WL: Physiology and electrical activity of uterine contractions. Semin Cell Dev Biol. 2007, 18 (3): 289295. 10.1016/j.semcdb.2007.05.004.View ArticlePubMedPubMed CentralGoogle Scholar
 La Rosa PS, Eswaran H, Lowery C, Preissl H, Nehorai A: Forward modeling of uterine EMG and MMG contractions. IFMBE Proceedings 11th World Congress on Medical Physics and Biomedical Engineering, Volume 25/4. 2009, Edited by Dssel O, Schlegel WC. Munich: Springer Berlin Heidelberg, 16001603.Google Scholar
 Tung L: A bidomain Model for Describing Ischemic Myocardial DC Potentials (Ph. D. thesis). 1978, Cambridge, Massachusets: M.I.T, Electrical Engineering and Computer Science DepartmentGoogle Scholar
 Miller W, Geselowitz D: Simulation studies of the electrocardiogram, I. The normal heart. Circulation Res. 1978, 43 (2): 301315. 10.1161/01.RES.43.2.301.View ArticlePubMedGoogle Scholar
 Peters MJ, Stinstra JG, Hendriks M: Estimation of the electrical conductivity of human tissue. Electromagnetics. 2001, 21 (7/8): 545557.Google Scholar
 Sadiku M: Elements of electromagnetics. 2001, New YorkOxford: Oxford university pressGoogle Scholar
 Plonsey R, Heppner DB: Considerations of quasistationarity in electrophysiological systems. Bull Math Biophys. 1967, 29 (4): 657664. 10.1007/BF02476917.View ArticlePubMedGoogle Scholar
 Malmivuo J, Plonsey R: Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields. 1995, USA: Oxford University PressView ArticleGoogle Scholar
 Rihana S, Marque C, Lefrancois E: A two dimension model of the uterine electrical wave propagation. Engineering in Medicine and Biology Society, 2007. EMBS 2007. 29th Annual International Conference of the IEEE. 2007, IEEE, 11091112.View ArticleGoogle Scholar
 Wikland M, Lindblom B: Relationship between electrical and mechanical activity of the isolated termpregnant human myometrium. Eur J Obstet Gynecol Reprod Biol. 1985, 20 (6): 337346. 10.1016/00282243(85)900577.View ArticlePubMedGoogle Scholar
 Kawarabayashi T, Kishikawa T, Sugimori H: Effect of oxytocin on spontaneous electrical and mechanical activities in pregnant human myometrium. Am J Obstet Gynecol. 1986, 155 (3): 671676.View ArticlePubMedGoogle Scholar
 Parkington HC, Tonta MA, Brennecke SP, Coleman HA: Contractile activity, membrane potential, and cytoplasmic calcium in human uterine smooth muscle in the third trimester of pregnancy and during labor. Am J Obstet Gynecol. 1999, 181 (6): 14451451. 10.1016/S00029378(99)70390X.View ArticlePubMedGoogle Scholar
 Shmygol A, BruMercier G, Gullam J, Thornton S: Control of uterine Ca2+ by membrane voltage. Ann New York Acad Sci. 2007, 1101: 97109. 10.1196/annals.1389.031.View ArticleGoogle Scholar
 FitzHugh R: Impulses and physiological states in theoretical models of nerve membrane. Biophys J. 1961, 1: 445466. 10.1016/S00063495(61)869026.View ArticlePubMedPubMed CentralGoogle Scholar
 Nagumo J, Arimoto S, Yoshizawa S: An active pulse transmission line simulating nerve axon. Proc. IRE, Volume 50. 1962Google Scholar
 FitzHugh R: Mathematical models of excitation and propagation in nerve. H P Schwan, ed Biol Eng. 1969, Chapter 1: 185.Google Scholar
 Franzone PC, Guerri L, Tentoni S: Mathematical modeling of the excitation process in myocardial tissue: influence of fiber rotation on wavefront propagation and potential field. Math Bioscienses. 1990, 101: 155235. 10.1016/00255564(90)90020Y.View ArticleGoogle Scholar
 Peters MJ, Stinstra JG, Uzunbajakau S, Srinivasan N: Fetal Magnetocardiography. Adv Electromagn Fields Living Syst. 2005, 4: 140. 10.1007/0387240241_1.View ArticleGoogle Scholar
 Weber R, Dickstein F: On the influence of a volume conductor on the orientation of currents in a thin cardiac issue. Lecture Notes Comput Sci. 2003, 2674: 111121. 10.1007/3540448837_12.View ArticleGoogle Scholar
 Wolfs G, Leeuwen M: Electromyographic observations on the human uterus during labour. Acta obstetricia et gynecologica scandinavica. 1979, 58 (S90): 161. 10.3109/00016347909156375.View ArticleGoogle Scholar
 The prepublication history for this paper can be accessed here:http://www.biomedcentral.com/17566649/12/4/prepub
Prepublication history
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.