Estimating uterine source current during contractions using magnetomyography measurements

Understanding the uterine source of the electrophysiological activity of contractions during pregnancy is of scientific interest and potential clinical applications. In this work, we propose a method to estimate uterine source currents from magnetomyography (MMG) temporal course measurements on the abdominal surface. In particular, we develop a linear forward model, based on the quasistatic Maxwell’s equations and a realistic four-compartment volume conductor, relating the magnetic fields to the source currents on the uterine surface through a lead-field matrix. To compute the lead-field matrix, we use a finite element method that considers the anisotropic property of the myometrium. We estimate the source currents by minimizing a constrained least-squares problem to solve the non-uniqueness issue of the inverse problem. Because we lack the ground truth of the source current, we propose to predict the intrauterine pressure from our estimated source currents by using an absolute-value-based method and compare the result with real abdominal deflection recorded during contractile activity. We test the feasibility of the lead-field matrix by displaying the lead fields that are generated by putative source currents at different locations in the myometrium: cervix and fundus, left and right, front and back. We then illustrate our method by using three synthetic MMG data sets, which are generated using our previously developed multiscale model of uterine contractions, and three real MMG data sets, one of which has simultaneous real abdominal deflection measurements. The numerical results demonstrate the ability of our method to capture the local contractile activity of human uterus during pregnancy. Moreover, the predicted intrauterine pressure is in fair agreement with the real abdominal deflection with respect to the timing of uterine contractions.


Introduction
It is an inverse problem to estimate the underlying source currents, such as location and time courses, from electromagnetic measurements of uterine contractions. Solving this inverse a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Materials and methods
In this section, we describe the collection of real and synthetic MMG data, and discuss the source current distribution of uterine contractions during pregnancy. We introduce a forward model for the magnetic field based on a lead-field matrix that is constructed on a realistic fourcompartment volume conductor, and provide the estimation of the underlying source currents and the corresponding intrauterine pressure.

Clinical site and MMG data
Three real MMG data sets were used for the estimation of uterine source current. These data sets were collected from two pregnant women at the University of Arkansas for Medical Sciences (UAMS), after the study protocol was explained and written consents to perform the study were obtained. The protocol was approved by the UAMS Institutional Review Board. The SARA device (Fig 1) we used to non-invasively collect the abdominal MMG data is installed in a magnetically shielded room next to the labor and delivery unit in the UAMS, to reduce external magnetic fields which interfere with the biomagnetic field generated by human organs. This SARA system consists of 151 primary magnetic sensors spaced 3 cm apart (Fig 1a), arranged in a concave array, covering the maternal abdomen from the pubic symphysis to the uterine fundus, and laterally over a similar span (Fig 1b). Each sensor measures the magnetic fields at two magnetometers, one of which is close to abdomen and the other one is 8 cm away from the first one in a direction similar but not identical to be normal to the SARA surface (Fig 1c), and the sensor measurement is the difference between the measurements from these two magnetometers. The difference between surface normal and the actual sensor orientation is particularly notable in areas such as the lower central portion of the SARA device due to space constraints in placing sensors in its convex surface. The patient simply sits and leans forward slightly against the smooth surface of the array (Fig 1d), allowing the SQUID sensors to receive electrophysiological signals.
All of the MMG data sets were first recorded at 250 Hz and then downsampled to 32 Hz. For the preprocessing, we applied a 8th-order band-pass Butterworth filter of 0.1 − 1 Hz to attenuate interfering maternal and fetal cardiac signals. A 8th-order band-stop (notch) Butterworth filter of 0.25 − 0.35 Hz was also applied to suppress maternal breathing, which is a prominent signal around 0.33 Hz. Noisy sensors were then removed to avoid possible pollution for MMG measurements. Among these data sets, one set has simultaneous recordings of the abdominal deflection, which were collected using an air-filled bag that was placed between the maternal abdomen and the SARA system. During uterine contractions, the pressure on the airbag induced by the abdominal shape change was transmitted via a tube to a pressure sensor that was connected to a standard fetal monitor which was located outside the shielded room. The output of the monitor was digitized and synchronized with the MMG signals. This simultaneous recording was performed as a proof of concept study and is difficult for routine application since noise artifacts could be introduced in the MMG data due to application of an external device.
Synthetic MMG data sets were employed to test our estimation approach. The synthetic MMG data sets were generated using our realistic multiscale electromagnetic model which was proposed in [9]. In this model, the volume conductor was exactly the same as the one in this inverse estimation work, and a sensor model was used to replicate the true SARA sensor positions and sensing directions as illustrated in Fig 1c. In particular, we represent the volume conductor (Fig 2) as four compartments (from the inner layer to the outer layer: fetus, amniotic fluid, uterus, and abdomen) with electrically conductive boundaries between compartments (Fig 2c). The geometry of an anatomically realistic uterus (Fig 2a) is based on the magnetic resonance images (MRI) of a pregnant, near-term woman and a uterine mesh is adopted from the FEMONUM project [10]. Considering the anisotropic nature and fiber variations of the myometrium, we divide the entire uterus into 25 contiguous regions. We set the region centers via random sampling from the finite-element mesh of uterus and then divide the regions by resampling any point that lies less than 4 cm from its nearest center. For each region, the fiber orientation, with respect to uterine surface tangential vector, is sampled from a normal distribution N $ ð0; p=4Þ (see red arrows in Fig 2a for detailed fiber orientations). Assuming the cylindrical symmetry of fibers, the longitudinal and transversal conductivity values for the myometrium are 0.68 S/m and 0.22 S/m, respectively. We also assume that the abdomen (Fig 2b) deforms to follow the shape of the SARA device when patients lean against it to take the MMG measurement. The corresponding conductivity values for the abdomen, amniotic fluid, and fetus are 0.2 S/m, 1.74 S/m, and 0.5 S/m, respectively (see more details in our previous work [9]).
In the uterus, the myometrial cells either can spontaneously generate their own impulses, pacemaker cells, or can be excited by the action potentials propagating from the neighboring cells, pacefollower cells. However, it is still unclear whether there exists a specific pacemaker mechanism or a specific pacemaker area [11]. Except for some observations of contractile activity originating from specialized cells [12], most activities arise from any site throughout the myometrium [13,14]. Among these sites, fundus is one possibility [15][16][17], hence we set the initiation of electrical activity to occur at the fundus in the simulations. This is an example location of the initiation area where the activity is excited, which makes no difference in our following analysis on source current. The final synthetic MMG data was generated by adding white noise with 5 fV/ p Hz to the original synthetic time courses.

Generation of electromagnetic fields
The electromagnetic fields of uterine contractions can be derived from the quasistatic approximation of Maxwell's equations, since the frequency of the associated bioelectrical phenomena is typically below 1 kHz. Thus the time derivatives of the electromagnetic fields can be ignored as source terms. In the quasistatic approximation, where E and B are the electric and magnetic fields, respectively; 0 and μ 0 denote the permittivity and permeability of the free space, respectively; and ρ and J are the total charge density and current density, respectively. We divide the total current density, J(r), into two components: the volume current, J v (r) = σ (r)E(r), which is the result of the macroscopic electric field in the volume conductor, and the primary current, J p (r): where σ(r) is the macroscopic conductivity of the volume conductor. The primary current, J p (r), is related to the original biological activity, which is the source current density we consider in this work. From Eq (1), the electric field, E, can be represented as the negative gradient of a scalar electrical potential, V, as Therefore, the total current density in Eq (5) becomes From Eqs (3) and (7), we obtain that r Á (r × B) = 0 = μ 0 r Á J = μ 0 r Á (J p − σrV), hence which shows that the scalar electrical potential, V, can be solved analytically using finite element method (FEM) [18] given the source current density, J p , and proper boundary conditions.

Linear forward model
The forward problem in uterine contractions is to calculate the magnetic field, B(r), outside the abdomen generated by the source current density, J p (r), within the uterus. According to the Biot-Savart law, the magnetic field can be computed as where l = r − r 0 is the vector pointing from the source point r 0 to the observation point r with magnitude l = klk 2 . Here, the prime refers to quantities in the source region. Since l/l 3 = −r(1/l) = r 0 (1/l), Eq (9) becomes wheren 0 is the unit normal vector pointing outwards the source surface. For the total current density that approaches zero sufficiently fast when the source point, r 0 , goes to infinity, Eq (10) becomes With Eq (7), Since r × (σrV) = rσ × rV + σr × rV = rσ × rV, With rσ × rV = V(r × rσ) − r × (Vrσ) = −r × (Vrσ), the magnetic field, B, is represented as Based on the equivalence between Eqs (9) and (11), we obtain from Eq (14) that According to Eq (8), the electrical potential, V, in the above equation can be computed from the source current density, J p , and is linearly related to J p . Therefore, based on Eq (15), the magnetic field, B, is linearly related to the source current density, J p [19]. Thus, there is a lead field, Lðr; r 0 Þ, relating the magnetic field measurement, B, at r to the source current, J p , at r 0 , satisfying If the source current, J p , is a current dipole with moment q = qd q in location r q , i.e., J p (r) = qδ(r − r q ), the magnetic field is given by where q = qd q denotes a current dipole with magnitude q pointing at direction d q .
For simplicity, the orientation of the current dipole is assumed to be perpendicular to the surface. When taking measurements using the SARA device, we cannot obtain all the components of the magnetic field; instead, we get a specific sensor-oriented component. Under these conditions, Eq (17) becomes where b denotes the sensor-oriented magnetic field measurement that is generated by a normal current dipole with magnitude q. Note that according to Eq (18), the lead field Lðr; r 0 Þ is exactly the same as the magnetic field measurement b at r if we apply a unit normal current dipole with q = 1 at r 0 . If the volume conductor is spherically symmetric and piecewise homogeneous, the lead field for the normal component of the magnetic field has a closed form [19]. However, it is a great challenge to obtain an explicit expression for the lead field in a complex volume conductor. We choose to solve this problem using FEM, which is a powerful tool for numerically solving for the lead field since it can deal with the anisotropy and realistic geometry in our volume conductor. In this work, considering the random initiation area [13,14], we adopt the distributed source current instead of a small number of current dipoles. In this case, a much larger number (usually greater than 5,000) of current dipoles are distributed over the whole myometrium. Since the myometrium is a thin layer [20], it is feasible for us to assume that the source current is limited to the external uterine surface. We divide the external uterine surface into small elements and introduce a current dipole at every vertex of the elements. Therefore, the lead field Lðr i ; r 0 j Þ is the numerical solution b(r i ) of the sensor-oriented magnetic field at sensor location r i generated by a unit current dipole at r 0 j , i.e., where M is the number of SARA sensors and N denotes the number of current dipoles. The linear relationship between the current dipole amplitudes q t and the measurements b t at time t, therefore, is given by . ., M is the magnetic field measured by ith sensor at time t; 2. q t ðr 0 j Þ; j ¼ 1; 2; . . . ; N is the amplitude of the jth current dipole at time t; Ng is the lead-field matrix, with gðr i ; r 0 j Þ obtained using Eq (19).

Inverse estimation of source currents
The data collected by all sensors at time t can be expressed as for t = 1, 2, . . ., T, where e t is the measurement noise at time t. With the constructed lead-field matrix G, we are interested in estimating the current dipole amplitudes q t from the magnetic fields b t measured using the SARA device. In the MMG inverse estimation problem, the number of unknowns, N, is usually greater than 5,000, but the number of measurements, M, is about one hundred. Because of this ill-posed nature, the inverse problem has an identifiability issue, in that there is no unique mathematically correct solution for the problem. To resolve this issue, it is necessary to impose additional constrains on the current dipoles. Here, we are interested in the most significant source currents, hence we embed extra information as an ℓ 1 norm on the current distribution [21]. The resulting convex optimization problem is known as the Lasso problem [22]: where λ is a regularization parameter balancing the least-squares error and the ℓ 1 penalty. A small λ puts more emphasis on the least-squares error, whereas large values of λ emphasize the ℓ 1 penalty. Applying convex optimization algorithms [23], we can obtain the source current amplitudes after solving this Lasso problem.

Prediction of intrauterine pressure
During pregnancy, the uterine contractile activities appear in a form of an intrauterine pressure increase. We employ an absolute-value-based method [24] to predict the intrauterine pressure from our estimated source currents. To each estimated source current at each location, we first apply a 4th-order low-pass Butterworth filter with a cut-off frequency of 2 Hz and then downsample the source currents with a sampling frequency of 4 Hz to lower the computational complexity. After rectifying the source currents, we obtain the approximate energy by summing over all locations. We then smooth the energy by applying a 4th-order low-pass Butterworth filter with a cut-off frequency of 0.02 Hz and resample at the original sampling frequency, resulting in our predicted intrauterine pressure.

Results
In this section, we first illustrate the constructed lead-field matrix and then present numerical examples using both synthetic and real MMG data to illustrate our approach. For the construction of the lead-field matrix, we divided the external uterine surface into 12,412 vertices, i.e., N = 12,412. The sensor-oriented magnetic field is measured on the abdominal surface using the SARA device with 151 sensors, i.e., M = 151. The constructed lead-field matrix G is therefore a 151-by-12,412 matrix, corresponding to the SARA device with 151 sensors and 12,412 unit current dipoles on the external uterine surface.

Validation of the constructed lead-field matrix
We displayed the lead fields corresponding to particular source locations on the SARA device to test the feasibility of the constructed lead-field matrix G. In order to present the relationship between the lead field and source location and distance, the specific source locations were chosen to be pairs, such as cervix and fundus, left and right, front and back (from the front perspective). The corresponding results are illustrated in  right (0.106, −0.0004, 0.098) and the corresponding lead fields. We observe that the magneticfield patterns appear at the corresponding surrounding areas of the sources. The lead fields generated by the sources at the front (0.001, −0.075, 0.148) and back (−0.0001, −0.075, 0.014) of the uterus are presented in Fig 5. Note that the magnetic-field intensity generated by a source at the front is higher than that generated by the one at the back, which is in agreement with the inverse square nature of magnetic field with respect to the distance between source and observation points.

Estimation using synthetic MMG data
We utilized three synthetic MMG data sets (see details in Table 1 and Fig 6) to test our inverse estimation approach. Since fundus is one of the sites where uterine contractions are observed to arise [15][16][17], we excited the uterine activity at the upper left (from the rear perspective) of the uterus in the simulations. The parameters of the cell-level model in our multiscale model were set to generate plateau-type and bursting-type action potentials, which are the two predominant types in both a single uterine smooth muscle cell and isolated strips of myometrium [25,26]. In each simulation, the sampling frequency was 10 Hz with simulation length of 10 s. We obtained the synthetic MMG data b t , t = 1, 2, . . ., 100 as 100 151-by-1 vectors and the leadfield matrix G as a 151-by-12,412 matrix, i.e., M = 151, N = 12,412, and T = 100.
The initiation area of the first two synthetic MMG data sets is illustrated in Fig 7, in which the uterus is drawn in blue with the initiation area highlighted in green. The contour of the SARA device is sketched in black curves. The first data set represents the short-time oscillations of magnetic field (Fig 6a), that correspond to plateau-type action potential recruiting a small region in the upper left of the uterus (Fig 6c). In the second data set, local activity of bursting-type action potential (Fig 6f) was adopted to generate long-time oscillations of magnetic field (Fig 6d). Regarding the change of initiation areas during a single contraction or successive contractions [14], we shifted the initiation area from the location illustrated in In the literature, there are many ways to define the regularization parameter λ [27]. In general, the degree of the regularization should be consistent with the level of noise involved in the measurements. Therefore, the regularization parameter λ was set to be a little bit greater  than the noise in our estimation. The choice of the regularization parameter is beyond the scope of this work. Fig 9 shows the time courses of the estimated source currents for Data Set 1 in Table 1. Here, we were interested in the ones with higher intensity, hence set the threshold of the absolute value of amplitudes to be 2 × 10 −5 . Fig 10a illustrates snapshots of the estimated source currents on the uterine surface after thresholding at different time instants t = 0.5 s, 1.0 s, 2.0 s. The synthetic source currents generated using our multiscale model are presented as the Arrow Surface (red arrows on the uterine surface) in Fig 10b, in which the synthetic MMG on the uterine surface is also displayed. Red arrows in this figure reflect the direction of the source current. We observe that the synthetic source currents are distributed in a small local region in the upper left of the uterus at the beginning, then appear in a larger constrained neighboring region, and finally return to quiescence. Comparing Fig 10a with Fig 10b, we can see that the distribution area of the estimated source currents resembles that of the synthetic ones, although it is not exactly the same. Furthermore, while we did not consider the tangential component of the source current, the estimated source currents capture the emergence of local activities and the involvement of a larger excited area in the following contractile activities, which are in agreement with the tissue recruitment and contraction coordination via limited action potential propagation in local contractions.
The snapshots of the estimated source currents for Data Sets 2 and 3 in Table 1 are illustrated in Figs 11 and 12, respectively. We observe that similar results are obtained for both data sets, i.e., the occurrence of local activity and recruitment of neighboring area for the synthetic and estimated source currents are fairly well matched. Note from the results at different time instants, t = 2.8 s and t = 7.0 s, in Fig 12 that the estimated source currents reflect the

Estimation using real MMG data
We validated our inverse estimation method using three real data sets (see details in Table 2). Representative MMG signals after preprocessing that is described in subsection Clinical site and MMG data are illustrated in Fig 13. The MMG signals over the SARA sensors that covered one uterine contraction for 60 s are shown in Fig 13a. The signals from sensors MLE1 and MLF1 with high amplitudes are highlighted in red. We can see strong uterine activities in the lower left region of the abdomen. For the first MMG data set that lasts for 12 minutes, we removed MMG data from two noisy sensors and also removed two corresponding rows of the lead-field matrix, thus   Table 2). The temporal courses of the MMG data over 149 sensors are shown in Fig 14a, with amplitude in Tesla (T). The MMG data collected when the woman experienced no contractions is at the level of several 10 −13 s. Therefore, the regularization parameter λ was set to be 10 −12 , slightly greater than the noise, in this estimation. The estimated MMG and source currents are shown in parts b and c of Fig 14. As can be seen from Fig 14, the source-current patterns match well with the real MMG data, i.e., there are stronger source currents underlying For the other two data sets, we also removed the noisy sensors and the corresponding rows of the lead-field matrix before estimation (see Data Sets 2 and 3 in Table 2). In Figs 16 and 17, we illustrate the temporal courses and snapshots of the two MMG data sets and the estimated source currents, respectively. We observe that the results are similar to those of the first data set, i.e., the estimated source current density matches well with the local contraction in the lower region and upper left region of the uterus, respectively.

Intrauterine pressure prediction
Currently, intrauterine pressure catheter or tocodynamometer is widely used by physicians for assessing uterine contractions during pregnancy. To show the potential clinical implications of our source current estimation despite the lack of ground truth, we predicted intrauterine pressure for Data Set 1 in Table 2, which includes simultaneous measurement of abdominal deflection during a contraction as described earlier. Fig 18 shows the real abdominal deflection data in red and our predicted intrauterine pressure in blue. We observe that 83.33% of the predicted intrauterine pressure peaks display good predictive timing of uterine contractions when compared with the abdominal deflection peaks. Note that the intrauterine pressure predicted from the myometrial source currents is several seconds in advance of the measured abdominal deflection, which is in agreement with the fact that the uterine electrical activities induce the increase of intrauterine pressure. According to the estimated source currents, there exist local contractile activities in uterus. Based on the mechanotransduction mechanism proposed in [28], a local contraction slightly increases the intrauterine pressure, resulting in a global wall tension increase and the induction of more local contractions that generate high intrauterine pressure. One possible reason for the difference of phase shifts between the real and predicted contractions is that the predicted contractions are calculated based on the local activities while the real ones are for the global change. The phase shift is dependent on the complex tissue recruitment and contraction coordination following the local contraction. The predicted intrauterine pressure curves for Data Sets 2 and 3 in Table 2 are also presented in Fig 19.

Discussion
In this work, we tested our inverse estimation approach using synthetic MMG data sets which were generated using our previously developed multiscale model. The initiation areas were set at the fundus of the uterus (Fig 7) and the resulting electrophysiological activity further recruits more area of the uterus in local contractions (Figs 10 and 11). In real high-resolution recordings, there exist complex wave-propagation patterns, such as three or more wavefronts emerging at different positions of uterus and time instants and propagating in different directions [29]. According to our previous work [9], it is possible to obtain various magnetic-field patterns by changing the model configuration, such as initiation location, initiation time, and fiber orientations, which will affect the emerging position, emerging time, and propagating direction of the waves. MMG measurement is independent of any kind of reference, thus ensuring that the SARA device measures localized activity. The magnetic sensors of the SARA device are spaced 3 cm apart. Activities at different positions should be differentiable if the In our multiscale model, the coordination of uterine contractions is through the continuous propagation of action potential in smooth muscle cells. However, animal and human data indicates that action potentials propagate noncontiguously and over short distances [14,30,31]. We observe from the local contractions in real MMG data that the estimated source  Table 2.
https://doi.org/10.1371/journal.pone.0202184.g018 currents arise not only continuously, but also in areas that are not the neighbor of initiation area (Figs 15 and 17). One possible explanation for this is the mechanotransduction for longdistance signaling mechanism presented in [28]. In particular, one local contraction increases intrauterine pressure, then increases wall tension, inducing more local contractions to generate strong uterine contractile activity.
Uterine activities are generated by myometrial source currents [31]. It is common to model these currents using current dipoles [32]. According to the superposition principle, all complex sources can be approximated by multiple current dipoles. In this work, we are interested in the distribution of source currents and hence assume for simplicity that the current dipoles are perpendicular to the surface of the myometrium. Our numerical results (Figs 10-12, 15  and 17) show that the distribution of the source current estimation over the myometrium is in agreement with the synthetic source current and the MMG pattern measured using the SARA device, respectively. The snapshots of estimation results for synthetic MMG data (Figs 10-12) illustrate that our approach can track the emergence of local activity and the recruitment of larger area of source currents, although we do not consider the tangential component. Estimating the directions of source currents can be accommodated using the same framework. In this case, the inverse solution would be the amplitudes of the source currents in three orthogonal directions at each vertex, calculated using the lead-field matrix, in which each column corresponds to the lead field generated by a unit current dipole pointed in each direction.
Regarding the volume conductor, a set of concentric spheres with homogeneous and isotropic conductivities is the simplest model, in which case the radial component of magnetic field has a closed form with respect to a tangential current dipole [7,19]. A variant on the spherical model introduced in our previous work [33] is a set of spheres with the outer uterine layer shifted to the front of the abdomen. In [34], a more complex volume conductor represented by differently shaped ellipsoids is developed, which was constructed based on anatomic diagrams from Hunter's Anatomia Uteri Humani Gravidi [35]. In this project, we applied a more realistic volume conductor, proposed in our previous work [9], based on the MRI of a near-term woman and the abdomen that deforms to follow the SARA contour.
The spherical representation of the volume conductor geometry is a good approximation at the early stage of pregnancy, and its lead field can be expressed in a closed form. Our lead-field matrix, however, is constructed for a realistic geometry that is obtained after the acquisition of anatomical magnetic resonance images of a pregnant, near-term woman. In general, the computation of the lead field for this realistic geometry requires numerical solutions, for which we applied the finite element method, considering the anisotropic property of the myometrium. Numerical solution is computationally expensive and requires specifying conductivity values for each compartment. The conductivity values of the intracellular and extracellular domains, unlike those of the abdomen, amniotic fluid, and fetus [36], have not been reported. In this work, these values were obtained from the effective myometrium conductivity calculated by applying Archie's law [37]. Regardless, construction of a lead-field matrix by using precise experimentally-confirmed conductivity values remains desirable.
The performance of the inverse calculation is sensitive to the regularization parameter λ. In general, the degree of the regularization should be consistent with the level of noise in the measurement data. Accordingly, the choice of λ is often determined by popular methods such as the discrepancy principle, χ 2 test, L-curve, and generalized cross validation [27]. In this work, a fixed value of λ, determined according to the noise, was used for the estimation. Although we did not choose an optimal value, λ was set to ensure that the signal to noise ratios (SNRs) of the MMG data sets were in a consistent range in order to mitigate their potential influence. We postulate that optimizing the regularization parameter will further improve the performance of our estimation, the next natural research topic.
Source estimation has wide application in many different anatomical domains, such as the brain, heart, and uterus [5,[38][39][40][41]. The inverse problem is to estimate a large number of current dipoles from about one hundred measurements, which is mathematically ill-conditioned in the sense that various source configurations can produce the same magnetic field pattern. To solve it, it is necessary to impose additional constrains on the current dipoles. Among them, the most commonly used is the minimum-norm constraint, which imposes a weighted ℓ 2 norm on the source current distribution [19,[42][43][44]. A nonlinear smooth constraint is included using Bayesian methods [45], whose performance depends greatly on the choice of prior distributions. In spite of these methods, it is quite difficult to validate the estimation accuracy. However, our first attempt at source estimation for uterine contractions, despite its limitations, is promising. Future research on developing a method to solve this inverse estimation is needed to achieve good estimation performance.

Conclusions
We proposed a method to estimate the underlying source currents in the myometrium from noninvasive abdominal MMG measurements of uterine contractile activities during pregnancy. Our method incorporates electrophysiological and anatomical knowledge of uterine contractions. We developed a forward model to describe the relationship between the abdominal magnetic field and myometrial source currents based on a lead-field matrix and used this model to compute the unknown source currents in the myometrium. We introduced a realistic four-compartment geometry as the volume conductor model and a current dipole as the source model. We also applied the finite element method to construct the leadfield matrix. We obtained the estimation of underlying source currents in the myometrium by solving a constrained optimization problem. We also predicted the intrauterine pressure, which is clinically used for uterine contraction measurements, from the estimated source currents based on an absolute-value-based method. Finally, we displayed the lead fields that are generated by unit current dipoles at particular locations and illustrated our approach through numerical examples using both synthetic and real MMG data. We then estimated source currents and predicted the intrauterine pressure to show its clinical implications. Our work is potentially important as a starting point for helping characterize underlying activities of uterine contractions during pregnancy and future diagnosis of contractile dysfunction.