Inference of mechanical states of intestinal motor activity using hidden Markov models

Background Contractions and relaxations of the muscle layers within the digestive tract alter the external diameter and the internal pressures. These changes in diameter and pressure move digesting food and waste products. Defining these complex relationships is a fundamental step for neurogastroenterologists to be able define normal and abnormal gut motility. Results Utilising an in vitro technique that allows for the simultaneous recording of intraluminal pressure (manometry) and gut diameter (video) in an isolated section of rabbit colon, we developed a technique to help define the mechanical states of the muscle at any point in space and time during actual peristaltic movements. This was achieved by directly relating the changes in pressure to the changes in diameter along the length of the gut studied. For each individual measure of pressure or diameter, 3 dynamic state components were identified; increasing or decreasing changes or a stable period. Two additional static state components, fully contracted and fully distended, were defined for the diameter. Then qualitative mechanical states of the muscle activity were defined as combinations of these state components. A hidden Markov model was used to correlate adjacent-in-time samples, and the Viterbi algorithm was used to infer the most likely sequence of mechanical states based on the observed data. From this a spatiotemporal map of the mechanical states was produced, showing the regions of active contractions, active relaxations, or passive states along the length of the gut throughout the entire recording period. Conclusions The identification of mechanical muscles states based on gut diameter and intraluminal pressure was possible by modelling muscle activation with a hidden Markov model.


Background
The gastrointestinal tract extends from the mouth to the anus. Swallowed food travels along the length of this tube at an appropriate speed to allow for the breakdown and absorption of nutrients. The movement of content results from a complex series of muscular contractions and relaxations. These movements of the gut wall alter the pressure profiles within the gut [1] which in-turn cause the digesting contents to move in an oral or anal direction. Abnormalities in these motor patterns are associated with several prevalent and unpleasant disorders that cost health care systems billions of dollars per year [2]. The ability to define the relationship that exist between wall motion, intraluminal pressure and the flow of content is a fundamental step in understanding how altered motor patterns effect the transport of luminal content. However, accurately defining these relationships in the human gut in vivo is problematic; ethical constraints prevent detailed examination of real-time movement of the gut wall. To overcome this problem we developed an in vitro animal preparation that allowed us to record simultaneously, both intraluminal pressure (high-resolution manometry) and gut diameter (video) in real time, across varying length (15-90 cm) of intestine [3] (Figure 1a,b).
In the next stage of this research we published a theory based paper in which we developed a strategy based on simple principle of biomechanics to deduce the http://www.biomedcentral.com/1472-6793/13/14 Figure 1 An overview of data acquisition, and how mechanical muscle states relate to the data. (a) A schematic of the experimental setup. A section of rabbit colon is kept alive in an organ bath and simultaneous intraluminal pressure (fibre-optic catheter, sensors spaced at 10 mm) and diameter (camcorder) are recorded during peristaltic movements. (b) from these recordings a diameter (magenta) and pressure (green) spatiotemporal map is created. (c) at any location within the diameter-pressure map we can plot the values of pressure (green) and diameter (magenta). In this example the values are taken from the 5s white line on the diameter-pressure map. (d) Plotting the relation between pressure and diameter in (c) is then displayed as an orbit plot. From the orbit plot the mechanical states of the muscle can be determined. In the example shown here 3 different mechanical states can been observed during the selected 5s period. The first mechanical state (blue-aqua section of the line) shows a decrease in diameter with no change in pressure (isotonic contraction; itc). In the next mechanical state (aqua-yellow) the pressure begins to increase with a corresponding decease in diameter (auxotonic contraction; atc). In the final mechanical state (yellow-red) the pressure drops with no further change in diameter (occluded isometric relaxation; oimr). mechanical state of the muscle (active contraction or relaxation, passive dilation, periods of quiescence, etc.) by calculating the relation between pressure and diameter at every point along the gut segment, and establishing where and when the muscle is actively contracting or relaxing [4] (Figure 1c,d). In that paper, space prevented a detailed description of the mathematical model involved in developing this process. In this work we address the details of the strategy that enabled us to confidently identify, for the first time, the mechanical states of the muscle during peristaltic contractions, and plot them as a functional spatiotemporal map.
A major factor in our work involved the use of hidden Markov models (HMM). The HMM have found extensive use in a variety of fields, stemming from their seminal application in speech recognition in the 1970s [5]. In gastroenterology, HMMs have been used to classify the location of a video capsule in the gastrointestinal tract [6]. In their application to hand written recognition, HMMs were utilized by observing lines and curves drawn on a 2D plane and inferring the most probable characters they represent [7]. In developing our intestinal model of mechanical states, one of the major challenges we faced was how to identify the beginning and the end of the periods of the various different states of the muscle. However, analogous to hand written recognition, our diameter-pressure plots can be observed as lines in the 2D plane. Therefore the HMM was adapted to infer the most probable muscles states at any given point in time along the length of the studies gut segment. A description of the steps involved in this process are detailed here.

Identification of muscle states
The circular muscle of the gut is arranged as continuous rings of smooth muscle forming the inner muscle layer of the tubular structure. The circular muscle layer is thicker than the external longitudinal muscle layer and thus is mechanically more powerful. When it contracts, propulsion of the luminal content can occur. This prolusion is made more efficient if the adjacent regions of the gut (proximal or distal) concurrently relax, thus effectively http://www.biomedcentral.com/1472-6793/13/14 allowing the bolus to move into a region of lower pressure [8]. These regions of contraction and relaxation can be detailed by simultaneously recording the diameter of the gut and the intraluminal pressure ( Figure 1a). From the resultant spatiotemporal map of diameter and pressure ( Figure 1b) we can plot, at any point in time, pressure and diameter curves (Figure 1c). From these we can then begin to express the mechanical relation between pressure and diameter as a pressure-diameter orbit plot (Figure 1d). In publishing the theory behind this concept [4] we also established that twelve distinct mechanical muscle states can be predicted ( Figure 2). The mechanical states are summarized below; Occluded isometric contraction (oimc): active; occurs with an increase in pressure, with no change in diameter, when the gut is at its minimum diameter. Occluded isometric relaxation (oimr): passive; occurs with a decrease in pressure and no change in diameter, when the gut is at its minimum diameter. Distended isometric pressure increase (dipi): passive; occurs with an increase in pressure, no change in diameter, when the gut is at its maximum diameter. Distended isometric pressure decrease (dipd): passive; occurs with a decrease in pressure, no change in diameter, and a maximum diameter.
Isotonic contraction (itc): active; occurs with a decrease in diameter and no change in pressure.
Isotonic relaxation (itr): active; occurs with an increase in diameter and no change in pressure. Auxotonic contraction (atc): active: occurs with a decrease in diameter and an increase in pressure. Auxotonic relaxation (atr): active: occurs with an increase in diameter and a decrease in pressure.
Passive shortening (ps): passive; occurs with a decrease in diameter and pressure.
Passive dilation (pd): passive; occurs with an increase in diameter and pressure.
Occluded quiescence (oq): passive; occurs with no change in diameter nor pressure, when the gut is at its minimum diameter. Distended quiescence (dq): passive; occurs with no change in diameter nor pressure, when the gut is at a non-minimum diameter.
Having established that 12 mechanical states exist we then need to develop a means to determine when and where each state existed along the gut in an automated fashion. In addition we had to detect both the changes in the mechanical states (transitions of direction of the linear segments of the orbits) and the permanence in that state (the beginning and end of the linear segment).

Methods
The intraluminal pressures were recorded by a manometry catheter with pressure sensors spaced at 10 mm intervals and each of these was considered independent for modelling purposes. At each sensor location, diameter and pressure values were recorded continuously.

Data acquisition
A detailed description of the techniques used to collate the data has been published elsewhere [4]. Here we will provide a brief overview.
Rabbits were euthanized humanely by intravenous injection of pentobarbitone sodium (0.5 ml kg −1 ) in accordance with approval by the Animal Welfare Committee of Flinders University.
The gut diameter and internal pressure were recorded using two different methods, at different sampling rates, spacial resolutions, and offsets in time and space. Resampling and spatiotemporal alignment was required to coincide values of pressure and diameter for use as a single observation vector for each node in the Markov chain. Each spatial location was assigned its own Markov chain.

Figure 2
Conceptual directions of observed data and their relation to mechanical muscle states. The left and middle graphs correspond to states where myogenic, or neurogenic excitatory or inhibitory mechanisms are active. The graph on the right corresponds to passive states where no muscle activation is required to exhibit that behaviour. http://www.biomedcentral.com/1472-6793/13/14 A spatiotemporal map of colonic diameter was based on techniques developed in our lab [9]. Briefly the spatiotemporal maps were obtained by recording a top-down video of the colon suspended in a bath of Krebs solution, such that the length of the colon appeared horizontally in the video, and the number of pixels that the colon spans for each vertical pixel line in each frame was counted. Using a reference ruler visible in the recorded video, the colonic diameter in millimetres was obtained.
Pressures were recorded by 10 mm-spaced sensors at 10 Hz with a catheter inserted into the colon. Baseline drift was removed with iterated Gaussian minima smoothing [10].
Different spatiotemporal resolutions and offsets required resampling and alignment. A diameter map and a pressure map were combined by creating a grid of coordinates using the coarsest resolution in time and space from either map aligned with an adjustable spatiotemporal offset, and then binning the original maps into the grid. This resulted in a spatial and time resolution of 10 mm and 0.25 s respectively.
The alignment resulted in a single set of spatiotemporal coordinates that map to coincident diameter and pressure values, allowing diameters and pressures to be quickly compared for analysis without further interpolation. The temporal offset was determined by aligning events that were synchronously recorded by video as flashes from a light bulb (for diameter data) and embedded into the pressure recording as meta-data. The spatial offset was manually obtained by observing overlapping images of the two maps, while adjusting the offset to a value best representative of correct alignment that correlates with the video.
The entire data set consisted of 6 rabbits recorded with approximately 20 sensors over periods of 10 to 20 minutes each, resulting in just over 420,000 samples (almost 30 hours); orbits were created for each sensor.

Hidden Markov model
Given a sequence of observations of diameter and pressure at a single location in the colon, our objective was to infer the most-likely sequence of mechanical muscle states that could have resulted in those observations. Each mechanical muscle state was directly represented as a hidden state in the hidden Markov model. In our application, the dependence of a state on the previous state in the Markov chain was required as an implicit smoothing technique.

Observations
For each sensor location there are two values observed, the diameter d and pressure p, sampled at a frequency of 4 Hz. Sample i of the diameter and pressure is given by d i and p i , where i ∈ {1 . . . N}. The time derivatives of those values are denoted with an overdot and estimated using central differences aṡ where t = 0.25s for the 4 Hz sampling, andṗ defined in the same manner asḋ.
Observation i is defined as the vector will be written as x a:b for the sake of brevity.

States
We segment a sensor's recording and classify each observation into one of twelve discrete mechanical states of the set S = {oimc, oimr, dipi, dipd, itc, itr, atc, atr, ps, pd, oq, dq} The conceptual directions and positions of subsequences of observations on a diameter-pressure plot representing examples of each of the twelves states is depicted in Figure 2 and described in Figure 1a. The stateobservation model is based on developing a quantitative model for classifying the dynamics depicted in that figure.

State-observation model
The state-observation model defines the conditional probabilities of observations and states. The goal is to arrive at a formula describing the distribution of possible observations, under the assumption that they were produced while the muscle was in any given mechanical state. The sample subscript i is omitted for brevity in this section, since the state-observation model is independent of the sample number. The states in S can be considered a vector of state components (c d , cḋ, cṗ), where cḋ ∈ {quietḋ, posḋ, negḋ} cṗ ∈ {quietṗ, posṗ, negṗ} Given a velocity observation v ∈ {ḋ,ṗ}, the probability of the corresponding state component being quiet is modelled on the basis of a normal distribution located at 0 with width free parameters σḋ and σṗ. The normal distribution is the most common way to model noise and errors, justified by the central limit theorem. Since very low speeds are considered quiescent, slow tonic muscle activity will be absorbed into passive states, with only phasic contractions exposed as active states. The probability of the state component being pos or neg is modelled on the basis of normal cumulative distribution functions with the same parameters. Probability mass functions are normalised to ensure See (10)-(12) for the formulae, and Figure 3 (a) for a visual example.
Given an observation of diameter d, the probability of the corresponding state component being occ or dis is modelled on the basis of a normal cumulative distribution function and its complement, located at μ d + min d with width σ d , where min d is the value of the minimum recorded diameter. The value d − min d is referred to as "dilation", and μ d represents the amount of dilation that separates occluded and distended states. The probability of the state being any given d is a constant. The probability mass functions are normalised to ensure See (13)-(15) for the formulae, and Figure 3 (b) for a visual example. The probability mass functions (PMFs) f (c o |o) are given by where is the cumulative standard normal distribution function, v is an alias for eitherḋ orṗ, and Z v and Z d are normalisation constants ensuring the probabilities sum to 1 for any given observation and all corresponding states. The probability density function (PDF) of the model emitting a particular observation element o ∈ {d,ḋ,ṗ} given a corresponding state component c o , is known as the emission PDF, which can be obtained from Bayes' theorem The emission PDF is used to determine the distribution of observations that can be made assuming an underlying muscle mechanical state. Since only the relative sizes of emission PDFs are required for calculating the optimal state sequence (rather than the actual probability values themselves), the observation prior density f (o) can be removed which is independent of state and so does not change which state is most likely to have produced a given observation. In a similar way, if we approximate the state component prior distribution f (c o ) with a constant, then it can also be removed resulting in the following practical approximation Note that under this approximation, f (o|c o ) no longer represents a PDF, that is, f (o|c o )do = 1. This approximation isn't necessary for performance, since finding the priors numerically is straight forward with Gaussian kernel density estimation of f (o) and the .09 for our data. Rather it is our desire to keep the emission PDFs free of any other data-based attributes besides the few parameters that were chosen by hand. This allows for the consistent interpretation of results by having a homogeneous muscle mechanics model which is comparable among different preparations and recordings.
Assuming independence of the variables d,ḋ, andṗ, the joint emission PDF can be factorized into with implicit substitutions of (10)-(12) and (13)-(15) to be made in (18) due the approximation (17). The independence assumption is a model simplification rather than an observed independence in the data, and allows for the Figure 4 A matrix of state transition occurrence counts. The vertical and horizontal axes represent each of the 12 states, such that a transition moving from one state (vertical axis) to another state (horizontal axis) is represented as a square shaded from black (0 occurrences of that transition) to white (many occurrences of that transition). Transitions that never occurred are marked with a red circle, and transitions that occurred only very rarely (between 1 and 4 times for the entire data set) are marked with a yellow circle. http://www.biomedcentral.com/1472-6793/13/14 factorization of the joint emission PDF into component emission PDFs. Making such an independence assumption, where none exists in reality, is common practice in machine learning for simplifying models and algorithms used to make inference. In practice, inference algorithms are often robust despite such assumptions, and we observe that inferences obtained by our technique remain consistent with expected states.

Optimal state sequence
Given an observation sequence o 1:N , our goal is to find the optimal sequence of states s 1:N that explains the observations. The joint PDF of the observation and state sequences is given by a hidden Markov model [5], which subsumes the following factorisation To keep the number of free parameters small, we choose the state transition probabilities to be independent of state for non self-transitions, given by where γ ∈ (0, 1) is a free parameter that represents a penalty (when γ > 1 |S| ) for changing states between samples. Such a penalty allows the inferred state sequence to be desensitised to noisy observations, resulting in practical smoothing of the sequence. It is important to include this implicit smoothing so that linear segments of orbits containing undulations are not subdivided, which may result in incorrect state inferences based on arbitrary angles of small subsegments of the undulated linear segment. A value of γ = 1 |S| would result in modelling the states in adjacent samples as independent, effectively eliminating the links in the Markov chain.
Due to the difficulty in quantifying and translating physically realistic dynamics into transition probabilities under our model, state-dependent transition probabilities were not used. However, by plotting all of the transition occurrences from our complete data set (over 420,000 samples) we have been able to demonstrate that qualitatively unlikely transitions were either not observed or rarely observed in the inference results ( Figure 4). For example we did not have a single example of a mechanical state moving from occluded quiescence to a distended isometric pressure increase. Thus we are confident that state-independent transition probabilities are sufficient.
The probability of the initial state is modelled with the discrete uniform distribution The solution to the most likely sequence of states for a given sequence of observations is obtained by performing the following maximisation arg max , and the solution can be found by applying the Viterbi algorithm. The Viterbi algorithm works by calculating the most likely sequence of states based on a sequence of discrete observations. To handle continuous observations we used the value of the emission PDF at each observation as a substitute for the probability of making that observation.

Parameters
Through an iterative trial-and-error approach we manually selected and tuned parameter values, checking selected areas of data where we could verify whether the inference resulted in states consistent with existing domain knowledge of the expected muscular activity.
For example, when the gut contracts and squeezes the catheter we expect the pressure to rise. For this approach, data was drawn from a total of 6 rabbits, each with 10 to 20 minutes of recording. The values for parameters are given in the middle column of Table 1.
To test the classification sensitivity of different parameter values to the ones selected manually, parameters were varied one at a time through a range of values (right column in Table 1). One reference classification, considered as ground-truth, was performed with the manually selected parameter values (middle column in Table 1). The reference classification was compared to many comparison classifications with the varied parameter values. A ratio of the number of incorrect active-state In (e) the gut is color coded with the colors from the multiple mechanical states shown in (b). In (f) the simplified mechanical states from (c) are shown. The isotonic contraction at 190 mm (orange) and auxotonic contraction at 180 mm (yellow) push the bolus to the right, while the isometric contraction at 150-170 mm (red) ensures no back-flow. The isotonic relaxation at 210-240 mm (light-blue) help to expand the colon to accommodate the incoming bolus. The distended quiescence at 200 mm (light-green) depicts a momentary stillness where changes in diameter and pressure are trivial. The isometric relaxation between 80-140 mm (dark-blue) represents the recovery from earlier active contractions in that region that resulted in propelling the liquid bolus to the right. http://www.biomedcentral.com/1472-6793/13/14 classifications to the total number of states which were active in both reference and comparison classification over the entire spatiotemporal map was given by the equation  Figure 5) shows that the model is adequately robust to variations in the parameters, as the errors are qualitatively considered small. The test was performed on a 10-minute recording with 26 sensors of a single isolated rabbit colon which exhibited typical activity (for a combined 62400 samples in total).

Results and discussion
The optimal state sequence for each of the sensor locations (0-250 mm at 10 mm intervals) was independently inferred. When visualising, each state sequence was placed according to the location of the corresponding sensor in the colon, resulting in a 2D map where horizontal sequence strips are vertically stacked, shown in Figure 6 (a)-(c).
Parameters of the hidden Markov model were chosen manually, giving robust results with respect to variations in the parameters, shown in Figure 5. The γ parameter specifying transition probabilities should be adjusted to account for variations in time sampling frequency, which was not required here as our data consisted entirely of 4 Hz sampling.
The states inferred by our method coincided with hypothesised states based on manual observation of orbital plots and the expected mechanical function of the colon. This allowed qualitative analysis of colon dynamics as given in the description of Figure 6(d)-(f ). A more detailed description of these finding can be found in our previous paper [4].

Conclusions
The use of the hidden Markov model to discriminate mechanical states of the intestinal muscle in an isolated preparation of the rabbit colon has given experimental neurogastroenterologists a novel powerful tool to identify the active and passive states of the intestinal muscle.
The graphic representation of where the active contractions and relaxations occur in the intestine at any particular time ( Figure 6) will allow for the testing of many hypotheses currently proposed but not validated on the mechanisms responsible for the appropriate mixing and propulsive movements [4].
At the present the parameters described here have been shown to work with the rabbit distal colon. Whether or not the same criteria could be applied to different species with different sized diameters is still to be determined. It is likely that the parameter will need to be adjusted for each species (including human). We are currently setting up studies with different animal species to test this hypothesis.
Endnotes a [P] = 1 if P is true, 0 if P is false. b 1 A (x) = 1 if x ∈ A, and 0 otherwise.