Email updates

Keep up to date with the latest news and content from BMC Medical Physics and BioMed Central.

Open Access Research article

Multiscale forward electromagnetic model of uterine contractions during pregnancy

Patricio S La Rosa1, Hari Eswaran2, Hubert Preissl23 and Arye Nehorai4*

Author Affiliations

1 Department of Internal Medicine-General Medical Sciences Division, Washington University School of Medicine, Saint Louis, Missouri 63110, USA

2 Department of Obstetrics and Gynecology, University of Arkansas for Medical Sciences, Little Rock, Arkansas 72205, USA

3 MEG-Center, University of Tübingen, Tübingen, Germany

4 Department of Electrical and Systems Engineering, Washington University, Saint Louis, Missouri 63130, USA

For all author emails, please log on.

BMC Medical Physics 2012, 12:4  doi:10.1186/1756-6649-12-4

The electronic version of this article is the complete one and can be found online at: http://www.biomedcentral.com/1756-6649/12/4


Received:5 July 2011
Accepted:14 October 2012
Published:5 November 2012

© 2012 La Rosa et al.; licensee 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.

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 pre-term. 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 four-compartment volume conductor geometry. To model the current source density at the myometrium we use a bidomain approach. We consider a modified version of the Fitzhugh-Nagumo (FHN) equation for modeling ionic currents in each myocyte, assuming a plateau-type transmembrane potential, and we incorporate the anisotropic nature of the uterus by designing conductivity-tensor 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 adult uterus is a thick walled, hollow, muscular organ formed by three layers: the external serous perimetrium, the myometrium, and the inner mucous endometrium [1]. The non-pregnant uterus wall thickness is approximately 15 mm to 20 mm, and during pregnancy the uterine wall becomes thin with values at the 39th week of pregnancy ranging from 7.4 ±1.8 mm at the low anterior wall (lower segment) at the bladder interface, to 10.06 ± 1.9 mm at the posterior wall [3]. The myometrium is responsible for contractions, and it is formed by fasciculi which are comprised of sheet-like and cylindrical bundles of myocytes embedded in a connective tissue matrix [4]. The myocytes in a cylindrical bundle contract, thus shortening the smooth tissue and increasing uterus wall tension, hence increasing the intrauterine pressure. Figure 1 illustrates the microanatomy of the pregnant human myometrium.

thumbnailFigure 1. Diagram of microanatomy of pregnant human myometrium [4]. Red lines represent current flows.

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 3-dimensional myometrium for different initial conditions and stimuli. The authors use a reaction-difussion equation coupled with a Fizhugh-Nagumo 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 Ca2 + 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.

In this work, we propose a multiscale forward electromagnetic model of human myometrial contractions during pregnancy that jointly takes into account electrophysiological and anatomical knowledge at the cellular, tissue, and organ levels. Our model aims to help in the characterization of contractions and the prediction of labor using MMG [18] and EMG [19]. Here, we extend our partial results presented in [20]. Figure 2 illustrates the different levels considered in our modeling approach. In particular, our approach is twofold: first, we model the current source density at the myometrium, using models of myocyte electrophysiological activity and anisotropic conductivity, and second, we solve the forward electromagnetic problem, specifically, we compute the magnetic field and the action potential at the abdominal surface generated by the myometrial current-source density using Maxwell’s equations subject to a volume conductor geometry. To model the current source density at the myometrium we propose to apply a bidomain approach. The bidomain equations are a set of reaction-diffusion equations derived first for modeling the current sources of the myocardium as a function of the cardiac-myocyte transmembrane potential, and these equations proved to be a successful approach to study the functioning of the heart [21,22]. The diffusion part of the equations governs the spatial evolution of the transmembrane potential, and the reaction part is given by the local ionic current cell dynamics. Here we introduce a modified version of the FitzHugh-Nagumo (FHN) equation for modeling ionic currents in each myocyte. Though FHN does not consider explicitly the Ca2 + dynamics, the simplicity of the FHN model makes it an attractive candidate for modeling the propagation of depolarization waves in such large 2D and 3D simulations as the numerical examples presented in this work. We propose a general approach to design the conductivity tensor orientation for any uterine shape and estimate the conductivity tensor values in the extracellular and intracellular domains, using Archie’s law [23]. We illustrate our modeling approach using a spherical, four-compartment volume conductor geometry. As we elaborate later on, certain aspects of it stand as a fair approximation of the volume conductor geometry when performing MMG measurements.

thumbnailFigure 2. Illustration of the proposed modeling approach.

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 AT, and the identity matrix of size n×n is denoted In. The set <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M1','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M1">View MathML</a> denotes the vector space of symmetric n×nmatrices, and the spaces of nonnegative definite matrices and positive definite matrices are denoted by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M2','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M2">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M3','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M3">View MathML</a>, 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 four-compartment 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

Figure 3 illustrates the four-compartment volume conductor geometry for our problem, where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M4','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M4">View MathML</a> represents the abdominal cavity and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M5','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M5">View MathML</a> the boundary surface defined by the abdomen, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M6','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M6">View MathML</a> represents the myometrium, and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M7','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M7">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M8','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M8">View MathML</a> are its external and internal boundary surfaces, respectively. The volume denoted by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M9','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M9">View MathML</a> represents the space filled with the amniotic fluid that exists between the internal uterine wall <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M10','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M10">View MathML</a> and the boundary <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M11','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M11">View MathML</a> defined by the fetus volume <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M12','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M12">View MathML</a>. The vectors rand rindicate the positions of the observation point and source, respectively, with respect to the main axis of reference.

thumbnailFigure 3. Representation of the four-compartment volume conductor geometry and the forward electromagnetic problem of uterine contractions.

Extrauterine magnetic field and electrical potential

The electromagnetic properties of uterine contractions can be analized by solving a set of Maxwell’s equations [24,25] subject to boundary conditions given by the volume conductor geometry. In general, in a passive non-magnetic medium the total current density is the sum of the ohmic volume current and the polarization current, also known as the displacement current [24,25]. The ohmic current is the result of ions flowing in the medium and the displacement current is the result of a time varying electric field [24,25]. If the temporal variations of the electric field and the magnetic field, i.e. low frequencies, are small enough that the displacement current is very small with respect to the ohmic current, then it is valid to model the electromagnetic phenomena using the quasi-static approximation of the Maxwell Equations. Under the quasi-static approximation, changes in the sources that generate the electromagnetic field affect all field quantities instantaneously in the whole domain, and the total current density in the volume conductor geometry will be the sum of ohmic currents only. Biolectromagnetic fields vary slowly in time, with frequency components below 1KHz and with a spatial characteristic length scale several times much larger than the volume conductor of interest [26], i.e, much larger than the diameter of the uterus. Hence, changes in the bioelectric sources affect the bioelectromagnetic field in the volume conductor geometry instantaneously, which justifies the use of the quasi-static approximation of Maxwell’s equations [25,26]. Therefore, the extrauterine magnetic field B(r,t) at a position rand instant t is given as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M13','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M13">View MathML</a>

(1)

where μo is the permeability of free space and J(r,t) is the total current density (in A/m2). J(r,t)is given by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M14','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M14">View MathML</a>

(2)

where Js(r,t) is the uterine current density source and G(r)E(r,t)is the conduction current density (or return currents), as described by Ohm’s law, E(r,t) denoting the electric field established by Js(r,t)and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M15','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M15">View MathML</a> denoting the conductivity tensor defined by each compartment. Then, from the quasi-static conditions, ∇·J(r,t)=0, so ∇·G(r)E(r,t)=−∇·Js(r,t). Moreover, since ∇×E(r,t)=0, it follows that E(r,t)=−∇ϕ(r,t), where ϕ(r,t) denotes the potential. Thus, the equation that governs the relationship between electromyogram potentials and uterine current sources is

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M16','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M16">View MathML</a>

(3)

Solving the forward electromagnetic problem of uterine contractions requires the computation of B(r,t) and ϕ(r,t) at <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M17','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M17">View MathML</a> using Eqs.(1) and (3), assuming that Js(r,t)is known in <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M18','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M18">View MathML</a> and G(r) is known in all domains defined by the volume conductor geometry (see Figure 3).

The biological current sources Js(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 Js(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 current-source density model

In the myometrium, the intracellular and extracellular domains are both physically connected through membrane gates, and the intracellular domain is connected though gap junctions [1,19]. Therefore, we model the myometrium using the bidomain modeling approach. This approach represents the tissue (myometrium) as two interpenetrating extra-intracellular continuous domains with different conductivity values along and across the direction of the fiber [21,22], and it models the tissue using the generalized-passive cable equation. The bidomain modeling approach was originally derived for modeling the propagation of the transmembrane potential of the myocardium and proved to be a successful approach to study the functions of the heart [21,22]. Figure 4 shows a simplified illustration of the tissue and the bidomain approach, where ϕi(r,t) and ϕe(r,t) are the intracellular and interstitial potentials, respectively, and vm(r,t)=ϕi(r,t)−ϕe(r,t)is transmembrane potential. The conductivity tensors in the intracellular and extracellular domains are denoted by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M19','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M19">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M20','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M20">View MathML</a> (in S/m), and, using Ohm’s law, the current densities in each domain are given by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M21','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M21">View MathML</a>. The transmembrane volume current density in (A/m3) is denoted by jm(r,t) and is given by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M22','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M22">View MathML</a>

(4)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M23','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M23">View MathML</a>

(5)

where jion(r,t) is the ionic volume current density (in A/m3) of a myocyte, jstim(r,t) is the stimulus volume current density (in A/m3), cm is the membrane capacitance per unit area (in F/m2), and am is the surface-to-volume ratio of the membrane (in 1/m). To summarize the above description, the bidomain approach models the signal propagation in the tissue using the generalized passive cable equation model (Figure 4) which considers the current densities flowing through the extracellular domain given by Je(r,t), flowing through the intracellular domain given by Ji(r,t), and connecting both media given by jm(r,t).

thumbnailFigure 4. Illustration of the bidomain modeling approach.

Applying the conservation of charge to both domains, we obtain the following relationships:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M24','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M24">View MathML</a>

(6)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M25','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M25">View MathML</a>

(7)

Adding (6) and (7), we have ∇·(Ji(r,t) + Je(r,t))=0. Hence, the total current density in the myometrium is given by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M26','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M26">View MathML</a>

(8)

which can be expressed in terms of vm(r,t)and ϕe(r,t) as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M27','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M27">View MathML</a>

(9)

where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M28','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M28">View MathML</a> is the bulk myometrium conductivity tensor. Since spatial variations of vm(r,t)depend on the local establishment of a transmembrane current density, jm(r,t)≠0, we define the impressed current-density source as <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M29','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M29">View MathML</a>. Note that Js(r,t)exists only when the spatial gradient exists, i.e., only in a region where the myometrium is undergoing depolarization (excitation) or repolarization.

The total current at the myometrium J(r,t)depends on the spatio-temporal variations of vm(r,t)and ϕe(r,t), which are governed by the system of equations formed by Eqs. (5), (6) and (7). Using simple algebraic manipulations, the aforementioned system of equations can be written in terms of vm(r,t) and ϕe(r,t) only, obtaining the following equivalent expressions:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M30','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M30">View MathML</a>

(10)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M31','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M31">View MathML</a>

(11)

This set of reaction-diffusion 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 vm(r,t) and ϕe(r,t) depend on Jion(r,t), Jstim(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 traveling-wave solutions of these equations whose waveform depends on Jion(r,t) and whose initiation depends on Jstim(r,t). In what follows, we describe the models for both current densities, Jion(r,t)and Jstim(r,t), and for the conductivity tensors, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M32','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M32">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M33','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M33">View MathML</a>.

Ionic current model

The predominant type of transmembrane-potential waveforms measured in the pregnant human myometrium are spikes and plateau [1,28-31]. In this work, we focus on modeling the plateau-type transmembrane potential, as it has been more frequently observed [1,28-30]. In particular, we model Jion(r,t), using a variation of the FitzHugh-Nagumo (FHN) equations [32-34], as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M34','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M34">View MathML</a>

(12)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M35','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M35">View MathML</a>

(13)

where ε1, ε2, k, v1, v2, v3, δ, γ, and βare model constants, and w (in V) is a state variable of the model. The parameter ε1 (in Ωm2) 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 ε1is 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 v1, v2, v3 (in V), and k (in 1/V2) control the range of vm(r,t). Note that for a given set of k,v1, v2, and v3, 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 phase-space 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 Ca2 + dynamics, and, moreover, it assumes that changes in the intra- and extra-cellular 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

We also introduce a temporal-spatial model for Jstim, representing the stimulus due to pacemaker areas [1,19], as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M36','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M36">View MathML</a>

(14)

where hi(r,t) is a spatio-temporal function with range in [0,1], νi is the amplitude (in V), and Npis 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 pace-follower, 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 uterine-fiber 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 non-pregnant uterus by diffusion tensor magnetic resonance imaging (DTMRI). From the ex-vivo analysis of five non-pregnant 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.

Assume that the conductivity tensors are diagonal in a local coordinate system that is defined with respect to each myocyte and characterized by the unit vectors {e1,e2,e3}. In particular, Gi and Ge are diagonal matrices <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M37','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M37">View MathML</a> given by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M38','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M38">View MathML</a>

(15)

In order to take into account variable fiber orientation in the myometrium, we need to describe it in a global Cartesian coordinate system in which the local basis is defined at any point r as A=[a1(r) a2(r) a3(r)], where a3(r)is parallel to the main fiber axis. The representation of the tensors Giand Ge in terms of a global coordinate system is given by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M39','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M39">View MathML</a>

(16)

Assuming that the myocyte fiber conductivities in both domains are cylindrically symmetric, then σex=σey=σet, σix=σiy=σit, σiz=σil, and σez=σel. Therefore, the conductivity tensors can be expressed as follows [35]:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M40','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M40">View MathML</a>

(17)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M41','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M41">View MathML</a>

(18)

Hence, to construct the conductivity tensors as a function of r, it is enough to define the vector field a3(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

To design a3(r) at each point r, we represent the uterus as a hollow volume with uniform thickness, and we describe it as a union of mutually disjoint closed surfaces or layers. We use the implicit definition of a surface, namely, the set of points r satisfying f(r)=0. Then, at each point r, we define a set of local orthonormal coordinate axes given by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M42','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M42">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M43','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M43">View MathML</a>, where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M44','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M44">View MathML</a> is the normal vector to the layer containing r, and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M45','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M45">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M46','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M46">View MathML</a> are mutually orthogonal vectors that belong to the tangent plane of the respective layer at point r. We define <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M47','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M47">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M48','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M48">View MathML</a> using the curve of symmetry of the uterine inner-circular layer as a reference [5]. This curve goes from the fundus to the cervix, and it coincides with the long axis of the non-pregnant uterus. We denote C as the curve of symmetry using the following parametric representation as a function of a single parameter t:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M49','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M49">View MathML</a>

(19)

where rC(t)is a point defined with respect to the global coordinate system, and rC(t1) and rC(t2) are the extreme points of the curve. For example, let rC(t)=(x(t),y(t),z(t)) with respect to the Cartesian system. Define <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M50','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M50">View MathML</a>, the unitary vector field with direction given by the tangent vector of rC(t) at t0, where rC(t0) is the closest point to r such that <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M51','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M51">View MathML</a>. Then, we define <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M52','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M52">View MathML</a> to be contained in the plane formed by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M53','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M53">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M54','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M54">View MathML</a> as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M55','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M55">View MathML</a>

(20)

subject to the following conditions:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M56','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M56">View MathML</a>

(21)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M57','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M57">View MathML</a>

(22)

Therefore, replacing (20) in (21) and (22) we obtain the following system of equations:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M58','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M58">View MathML</a>

(23)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M59','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M59">View MathML</a>

(24)

Solving for all r, such that <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M60','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M60">View MathML</a>, we obtain <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M61','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M61">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M62','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M62">View MathML</a> Then <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M63','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M63">View MathML</a>, since by definition it is mutually orthogonal to <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M64','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M64">View MathML</a>. Hence, given <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M65','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M65">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M66','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M66">View MathML</a> we define a3(r) as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M67','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M67">View MathML</a>

(25)

where α is the fiber orientation angle with respect to <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M68','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M68">View MathML</a>. In order to take into account complex fiber orientations αcan be modeled as a spatial function defined over the domain of interest. Given our uterine volume assumptions, the points rC(t1) and rC(t2) are the only points that satisfy the condition <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M69','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M69">View MathML</a>. Since at rC(t1) and rC(t2) we cannot define <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M70','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M70">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M71','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M71">View MathML</a> using the curve C as a global reference, we set a3(r)=0, defining a point of isotropic conductivity. In Figure 5 we represent the fiber orientation a3(r) with respect to the local coordinate axes given by <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M72','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M72">View MathML</a>.

thumbnailFigure 5. Simplified illustrations ofa3(r)with respect to the local coordinates axis given by<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M75','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M75">View MathML</a>,<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M76','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M76">View MathML</a>. The blue plane contains the vectors <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M77','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M77">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M78','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M78">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M79','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M79">View MathML</a>, and it is perpendicular to the gray plane formed by vectors <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M80','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M80">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M81','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M81">View MathML</a> and a3(r). The orange plane is the cross section of the uterus perpendicular to the vector <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M82','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M82">View MathML</a>. The gray curve is the curve of symmetry rC(t)with rC(t1)and rC(t2)extreme points of the curve.

If we have a uterine volume such that C is parallel to the z axis, then a3(r) can be written as a function of <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M83','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M83">View MathML</a> as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M84','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M84">View MathML</a>

(26)

where

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M85','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M85">View MathML</a>

(27)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M86','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M86">View MathML</a>

(28)

with ∇jthe j-th component of the gradient. In the case of a spherical myometrium, a3(r) is given as by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M87','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M87">View MathML</a>

(29)

where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M88','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M88">View MathML</a> Note that, for α=0, the main axis of the fibers runs vertically from the fundus to the cervix.

Estimating myocyte fiber conductivities

To the best of our knowledge, values of the intracellular and extracellular conductivity tensors have not been reported for the human myocyte, and therefore, these have to be estimated. To model the extracellular conductivity values σel and σet, we assume a grid-type distribution of myocytes in the myometrium and use an estimate of the extracellular conductivity the human myometrium obtained by applying Archie’s law [23]. We describe myocytes as long cylinders with diameter dcell and axis length lcell, such that dcelllcellAssuming that myocytes are uniformly arranged in a cubical grid with length lT=lcell + 2Δe and whose cross section has sides dT=dcell + 2Δe,then we have σeland σet as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M89','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M89">View MathML</a>

(30)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M90','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M90">View MathML</a>

(31)

where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M91','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M91">View MathML</a> is the conductivity of the extracellular medium in the myometrium. <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M92','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M92">View MathML</a> can be computed using the effective myometrium conductivity <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M93','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M93">View MathML</a>, available in the literature, and Archie’s law [23] as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M94','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M94">View MathML</a>

(32)

where p is the fraction of the volume occupied by the myocytes and collagenous fibers in the tissue, and m is the so-called cementation factor, which depends on the shape and orientation of the myocyte in the tissue.

To compute the intracellular conductivities σiland σit, we assume that the intracellular and extracellular domains have equal anisotropy ratios, that is,

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M95','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M95">View MathML</a>

(33)

and thus we need to compute ς. We obtain an analytical expression for ςusing reported values of the propagation speed of a transmembrane potential waveform traveling on isolated tissue strips from pregnant human myometrium at term [28]. In particular, replacing (33), (12), and (13) in the bidomain equations (10) and (11), and solving vm(r,t) for a traveling wave solution vm(ξ·r−ct), where ξ is a unitary vector pointing along the main axis of the myocyte and c is the speed of propagation, we obtain the following expression for ς:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M96','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M96">View MathML</a>

(34)

where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M97','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M97">View MathML</a>. Further, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M98','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M98">View MathML</a>, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M99','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M99">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M100','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M100">View MathML</a> are the roots of the following polynomial in vm:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M101','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M101">View MathML</a>

(35)

with vmrdenoting the resting transmembrane potential of the human myocyte. Note that in order to have ς≥0, ε1 must satisfy the following inequality:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M102','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M102">View MathML</a>

(36)

Monodomain approximation and boundary conditions

Assuming an equal anisotropy ratio, Eq. (33), simplifies the solution of the bidomain equations (10) and (11) by decoupling them as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M103','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M103">View MathML</a>

(37)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M104','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M104">View MathML</a>

(38)

The above simplification is also known as the monodomain approximation of the bidomain equations, which, under suitable boundary conditions, allows the computation of vm and, thus, Js, independent from ϕe.

To set up boundary conditions for computing electrical potentials, we need to take into account the volume conductor geometry (see Figure 3). In particular, we have two bidomain-monodomain interfaces: one between the myometrium <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M105','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M105">View MathML</a> and abdominal volume <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M106','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M106">View MathML</a> and one between the myometrium and the intrauterine cavity <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M107','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M107">View MathML</a>. Therefore, we have the following boundary conditions: (i) continuity of the interstitial potential ϕe at the perimetrium surface <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M108','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M108">View MathML</a> to the abdomen potential <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M109','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M109">View MathML</a>, (ii) flow of the normal component of Jthat crosses over from the uterus to the abdominal medium, (iii) no flow of the normal component of Jsto the abdominal medium, (iv) continuity of the interstitial potential ϕe at the endometrium surface <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M110','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M110">View MathML</a> to the intrauterine cavity potential <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M111','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M111">View MathML</a>, (v) flow of the normal component of Jthat crosses over from the uterus to the intrauterine cavity filled with amniotic fluid, (vi) no flow of the normal component of Js to the intrauterine cavity, (vii) no flow of the normal component of Jthat crosses over from the abdominal cavity to air, and (viii) either no flow or else flow of the normal component of J that crosses over from the intrauterine cavity, filled with amniotic fluid, to the fetus, depending on if it is covered with vernix caseosa (λ=0) or not (λ≠0) [36]. These boundary conditions are summarized as follows:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M112','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M112">View MathML</a>

(39)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M113','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M113">View MathML</a>

(40)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M114','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M114">View MathML</a>

(41)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M115','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M115">View MathML</a>

(42)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M116','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M116">View MathML</a>

(43)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M117','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M117">View MathML</a>

(44)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M118','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M118">View MathML</a>

(45)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M119','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M119">View MathML</a>

(46)

where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M120','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M120">View MathML</a> 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 plateau-type, and any myocyte can play the roll of pacemaker or follower by means of a stimulus current.

Numerical computation of vm(r,t), ϕ(r,t), and B(r,t)

The computations of vm(r,t), ϕ(r,t), and B(r,t)are given by the following procedure:

Step 1: Solve for vm(r,t)using Eqs. (37), (12), (13), and (14) subject to boundary conditions (41) and (44), and to initial conditions given by

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M121','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M121">View MathML</a>

(47)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M122','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M122">View MathML</a>

(48)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M123','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M123">View MathML</a>

(49)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M124','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M124">View MathML</a>

(50)

Step 2: Solve for ϕe(r,t)in <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M125','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M125">View MathML</a> and ϕ(r,t) in <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M126','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M126">View MathML</a> and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M127','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M127">View MathML</a> using the solution of vm(r,t), computed in Step 1, and the following expressions:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M128','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M128">View MathML</a>

(51)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M129','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M129">View MathML</a>

(52)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M130','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M130">View MathML</a>

(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 vm(r,t), ϕe(r,t), and ϕ(r,t), obtained in Steps 1 and 2, and considering the following boundary conditions:

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M131','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M131">View MathML</a>

(54)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M132','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M132">View MathML</a>

(55)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M133','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M133">View MathML</a>

(56)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M134','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M134">View MathML</a>

(57)

<a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M135','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M135">View MathML</a>

(58)

where <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M136','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M136">View MathML</a> 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 <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M137','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M137">View MathML</a>, and <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M138','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M138">View MathML</a> describe the continuity of the magnetic field between domains, and the boundary condition at <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M139','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M139">View MathML</a>, the external boundary of <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M140','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M140">View MathML</a>, 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 64-bit processors at 2.3 GHz, with 32 GB RAM.

Results

In the following, we illustrate our modeling approach by considering the electrophysiological characteristics of the myometrium and a spherical volume conductor geometry. In Figure 6, we illustrate a spherical, four-compartment volume conductor geometry used in our numerical example. We define a spherical myometrium with a 16 cm radius measured from the center to <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M141','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M141">View MathML</a> and assume the uterine wall has a uniform thickness of 1 cm. We also consider a spherical fetus with a 12 cm radius concentric to the myometrium and fully covered with vernix caseosa, i.e., λ=0in (46). The abdominal compartment is also spherical, with a 21 cm radius shifted −3 cm from the center of the myometrium in the x axis. This shift is to take into account that the uterus is closer to the abdominal surface than to the dorsal surface. We set the coordinate axis of reference at the center of the myometrium.

thumbnailFigure 6. Four-compartment volume conductor geometry used in the numerical examples. (a) View of z-x plane, and (b) z-y plane. Each compartment is assigned a different color. The myometrium has a non-uniform color to denote that its conductivity is anisotropic.

The conductivity values for each compartment are given in Table 1. In particular, to compute the extracellular myometrial conductivity tensors, we use the average values for the uterine myocyte dimensions at term based on data reported in [1,4,9] (see Table 2). We assume the average human myocyte shape to be a long cylinder with a small cross section; therefore, we use a cementation factor m=4/3 (see [23] for more details on the computation of this factor). The volume fraction p occupied by myocytes and collagenous fibers in the myometrium is set to 0.6. In order to consider an average myometrial fiber architecture that contains both circularly- and obliquely-oriented fibers we choose the fiber orientation angle α to be 45o. Figure 7 illustrates the global structure of the myometrial fiber orientation for this angle.

thumbnailFigure 7. Geometry and fiber orientation in spherical myometrium given byα=45°.

Table 1. Conductivity values of the volume conductor geometry

Table 2. Myocyte dimensions and Archie’s law parameters

We select the model parameters of the ionic current model using phase-space analysis, using as a reference the average plateau-type transmembrane potential recorded from isolated tissue strips of human myometrium at term [13,30]. In particular, the average resting potential, considering the results reported from the 37th week of pregnancy onwards, is approximately −56 mV. The plateau has an average depolarization of −27±1mV that terminates in 0.9±0.2minutes by an abrupt repolarization to the resting level [13,30]. Table 3 gives the parameter values used in the numerical example. Note that we compute the surface-area to volume-ratio am using the myocyte dimensions in Table 2.

Table 3. Ionic current model 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 jstim with Np=1, ν1=2 V, and h1(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 current-source densities and their electromagnetic fields at the level of the abdomen, we use a three-step 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 <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M146','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M146">View MathML</a> 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.

Figure 8 shows several snapshots of the FEM solution for one pacemaker on the fundus of a spherical myometrium, assuming anisotropy as illustrated in Figure 7. Figure 8 (a)-(c) illustrates the transmembrane potential and source current density distribution at the myometrium, Figure 8 (d)-(f) shows the electrical potential at the abdominal surface, and Figure 8 (g)-(i) illustrates the magnetic field density at the abdominal surface. The magnetic field measured at the abdominal surface, BMMG (the magnetomyogram field), is proportional to <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M147','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M147">View MathML</a> the projection of B onto the normal vector of the abdominal surface, <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M148','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M148">View MathML</a>

thumbnailFigure 8. FEM solution at time instants t=10 [s], 36 [s], 55 [s] for one pacemaker on the fundus of a spherical myometrium, assuming anisotropy. (a)-(c) transmembrane potential and source current density distribution at the myometrium, (d)-(f) electrical potential at the abdominal surface, and (g)-(i) magnetic field density at the abdominal surface.

In Figure 9 (a) we illustrate the temporal response of the FEM solutions for the transmembrane potential at different elevations over time. It can be seen that a stable traveling waveform has been established as the shape remains the same. Also, the maximum depolarization is −16 mV, the average potential in the plateau area is −25 mV, and the transmembrane potential duration, before hyperpolarization, is approximately 35 s, which are fair approximations with respect to the average recorded transmembrane potentials discussed in [13,30]. Note that our ionic current model introduces hyperpolarization, which constrains the excitability of the cell, and thus consecutive contractions can only take place until vmreaches 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 300s.

thumbnailFigure 9. (a) Temporal response of the FEM solutions for transmembrane potential at different elevations; (b) Percentage of contracting myometrial volume as a function of time.

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 Jsis rotated at a certain angle from the main direction of the propagating transmembrane potential vm, and it is the transversal component of this current, parallel to the x-y plane, which generates the magnetic field <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M149','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M149">View MathML</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 <a onClick="popup('http://www.biomedcentral.com/1756-6649/12/4/mathml/M150','MathML',630,470);return false;" target="_blank" href="http://www.biomedcentral.com/1756-6649/12/4/mathml/M150">View MathML</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, floor-mounted 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 population-average 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 plateau-type 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 vm 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 ε2value 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 preterm-labor.

Conclusion

We proposed a multiscale-forward 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 vm 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 FitzHugh-Nagumo equation for modeling the ionic currents in each cell. Though our ionic current model does not consider Ca2 + 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 pear-shaped 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

aSARA was built in collaboration with VSMMedTech Ltd., Canada and is installed at the University of Arkansas for Medical Sciences (UAMS) Hospital.

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

PSL contributed to the model design, implemented the model, and drafted the manuscript. HE and HP contributed to the model design and guided the model validation. AN coordinated the study, contributed to the model design, and helped in preparing the draft the manuscript. All authors read and approved the final manuscript.

Acknowledgements

The study was funded in part by NIH Grant NIBIB/1R01 EB007264-01A2, National Science Foundation, Grant No. CCF-0963742, and by Grant DFG BI 195-50.

References

  1. Chard T, Grudzinskas JG: The Uterus. Cambridge: University Press; 1995. OpenURL

  2. Lamont R: Looking to the future.

    BJOG: Int J Obstet Gynaecol 2003, 110(20):131-135. Publisher Full Text OpenURL

  3. Degani S, Leibovitz Z, Shapiro I, Gonen R, Ohel G: Myometrial thickness in pregnancy: longitudinal sonographic study.

    J Ultrasound Med 1998, 17(10):661-665. PubMed Abstract | Publisher Full Text OpenURL

  4. Young R, Hession RO: Three-dimensional structure of the smooth muscle in the term-pregnant human uterus.

    Obstet Gynecol 1999, 93:94-99. PubMed Abstract | Publisher Full Text OpenURL

  5. Weiss S, Jaermann T, Schmid P, Staempli P, Niederer P, Caduff R, Bajka M: Three-dimensional 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:84-90. OpenURL

  6. Young RC: Myocytes, myometrium, and uterine contractions.

    Ann New York Acad Sci 2007, 1101:72-84. Publisher Full Text OpenURL

  7. Andersen HF, Barclay ML: A computer model of uterine contractions based on discrete contractile elements.

    Obstet Gynecol 1995, 86:108-111. PubMed Abstract | Publisher Full Text OpenURL

  8. Barclay M, Andersen H, Simon C: Emergent behaviors in a deterministic model of the human uterus.

    Reprod Sci 2010, 17(10):948-954. PubMed Abstract | Publisher Full Text OpenURL

  9. Young R: A computer model of uterine contractions based on action potential propagation and intercellular calcium waves.

    Obstet Gynecol 1997, 89:604-608. PubMed Abstract | Publisher Full Text OpenURL

  10. Rabotti C, Mischi M, Beulen L, Oei G, Bergmans J: Modeling and identification of the electrohysterographic volume conductor by high-density electrodes.

    Biomed Eng, IEEE Transac 2010, 57(3):519-527. OpenURL

  11. 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):e11037. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  12. Aslanidi O, Atia J, Benson A, van den Berg H, Blanks A, Choi C, Gilbert S, Goryanin I, Hayes-Gill 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:183-192. PubMed Abstract | Publisher Full Text OpenURL

  13. Bursztyn L, Eytan O, Jaffa AJ, Elad D: Modeling myometrial smooth muscle contraction.

    Ann New York Acad Sci 2007, 1101:110-138. Publisher Full Text OpenURL

  14. Bursztyn L, Eytan O, Jaffa AJ, Elad D: Mathematical model of excitation-contraction in a uterine smooth muscle cell.

    Am J Physiol - Cell Physiol 2007, 292(5):C1816—C1829. PubMed Abstract | Publisher Full Text OpenURL

  15. Rihana S, Marque C: Modelling the electrical activity of a uterine cell, a mathematical model approach. In Proc. The 3rd European Medical and Biological Engineering Conference. Prague; 2005. OpenURL

  16. Rihana S, Terrien J, Germain G, Marque C: Electrophysiological model of the uterine electrical activity.

    Med Biol Eng Comput 2009, 47(6):665-675. PubMed Abstract | Publisher Full Text OpenURL

  17. 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):e18685. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  18. Eswaran H, Preissl H, Wilson JD, Murphy P, Robinson S, Lowery C: First magnetomyographic recordings of uterine activity with spatial-temporal information with 151-channel sensor array.

    Am J Obstet Gynecol 2002, 187:145-151. PubMed Abstract | Publisher Full Text OpenURL

  19. Garfield RE, Maner WL: Physiology and electrical activity of uterine contractions.

    Semin Cell Dev Biol 2007, 18(3):289-295. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  20. La Rosa PS, Eswaran H, Lowery C, Preissl H, Nehorai A: Forward modeling of uterine EMG and MMG contractions. In IFMBE Proceedings 11th World Congress on Medical Physics and Biomedical Engineering, Volume 25/4. Edited by Dssel O, Schlegel WC. Munich: Springer Berlin Heidelberg; 2009:1600-1603. OpenURL

  21. Tung L: A bi-domain Model for Describing Ischemic Myocardial D-C Potentials (Ph. D. thesis). Cambridge, Massachusets: M.I.T, Electrical Engineering and Computer Science Department; 1978. OpenURL

  22. Miller W, Geselowitz D: Simulation studies of the electrocardiogram, I. The normal heart.

    Circulation Res 1978, 43(2):301-315. PubMed Abstract | Publisher Full Text OpenURL

  23. Peters MJ, Stinstra JG, Hendriks M: Estimation of the electrical conductivity of human tissue.

    Electromagnetics 2001, 21(7/8):545-557. OpenURL

  24. Sadiku M: Elements of electromagnetics. New York-Oxford: Oxford university press; 2001. OpenURL

  25. Plonsey R, Heppner DB: Considerations of quasistationarity in electrophysiological systems.

    Bull Math Biophys 1967, 29(4):657-664. PubMed Abstract | Publisher Full Text OpenURL

  26. Malmivuo J, Plonsey R: Bioelectromagnetism: Principles and Applications of Bioelectric and Biomagnetic Fields. USA: Oxford University Press; 1995. OpenURL

  27. Rihana S, Marque C, Lefrancois E: A two dimension model of the uterine electrical wave propagation. In Engineering in Medicine and Biology Society, 2007. EMBS 2007. 29th Annual International Conference of the IEEE. IEEE; 2007:1109-1112. OpenURL

  28. Wikland M, Lindblom B: Relationship between electrical and mechanical activity of the isolated term-pregnant human myometrium.

    Eur J Obstet Gynecol Reprod Biol 1985, 20(6):337-346. PubMed Abstract | Publisher Full Text OpenURL

  29. 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):671-676. PubMed Abstract OpenURL

  30. 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):1445-1451. PubMed Abstract | Publisher Full Text OpenURL

  31. Shmygol A, Bru-Mercier G, Gullam J, Thornton S: Control of uterine Ca2+ by membrane voltage.

    Ann New York Acad Sci 2007, 1101:97-109. Publisher Full Text OpenURL

  32. FitzHugh R: Impulses and physiological states in theoretical models of nerve membrane.

    Biophys J 1961, 1:445-466. PubMed Abstract | Publisher Full Text | PubMed Central Full Text OpenURL

  33. Nagumo J, Arimoto S, Yoshizawa S: An active pulse transmission line simulating nerve axon.

    Proc. IRE, Volume 50 1962. OpenURL

  34. FitzHugh R: Mathematical models of excitation and propagation in nerve.

    H P Schwan, ed Biol Eng 1969, Chapter 1:1-85. OpenURL

  35. 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:155-235. Publisher Full Text OpenURL

  36. Peters MJ, Stinstra JG, Uzunbajakau S, Srinivasan N: Fetal Magnetocardiography.

    Adv Electromagn Fields Living Syst 2005, 4:1-40. Publisher Full Text OpenURL

  37. 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:111-121. Publisher Full Text OpenURL

  38. Wolfs G, Leeuwen M: Electromyographic observations on the human uterus during labour.

    Acta obstetricia et gynecologica scandinavica 1979, 58(S90):1-61. Publisher Full Text OpenURL

Pre-publication history

The pre-publication history for this paper can be accessed here:

http://www.biomedcentral.com/1756-6649/12/4/prepub