Human Duodenal MMC Phase 3 Motility Model Based on Manometry Readings
3/7/2018
1. Introduction
Migrating Motor Complex (MMC): In the fasted state, human gastrointestinal motility become organized as a cyclical motor pattern MMC consisting of: quiescent phase 1; phase 2 showing an increasing number of intermittent but irregular and rarely propulsive contractions; and phase 3 activity that is a group of the largest amplitude peristaltic waves occurring at their maximum frequency, and the entire group of contractions migrates distally over a long distance in an organized fashion. Ample experimental data support a coupling between interdigestive GI motility and secretion. Phase 1 secretion is characterized by minimal exocrine pancreatic and gastric secretion. During phase 2 secretion increases. Peak gastric acid secretion and bicarbonate secretion into the duodenum occur coinciding with the start of phase 3.
Therefore, an insight into phase 3, the decisvie factor of gastrointestinal motility and secreation, helps us better predict drug transit and dissolution in gastrointestinal tract and subsequent drug absorption profile.
This study is primarily aimed at building up a model of duodenal MMC phase 3 motility behavior in the fasted state based on manometry data via time series analysis.
2. Data Exploration
- The original manometry data was sampled every 0.1s for 7hrs. For the purpose of course project, the data in this study was subsetted by extracting one value out of every 1s during the phase 3 activity of first MMC cycle after determination at mid-duodenum of one subject. The summary info and plot along the time course are displayed below.
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -13.000 -4.000 1.000 1.413 6.000 30.000## [1] "SD = 7.48"- One may first want to check the potential seasonality since phase 3 is known for regular contractions relative to the other two phases. There is an obvious pattern occurring between every 6 lags and the period of the most prominant frequency in the periodogram is also 6s, indicating a mainstream period of 6s per peak in phase 3, which is in concert with the literature claim “10-12 contractions per minute”.
## [1] "period = 6"After differencing by 6 lags, the seasonal pattern disappear in the ACF plot and time plot looks more regularized than before, suggesting adding a seasonal component to the model. 
3. ARMA Model
- A time series data can be decomposed into three components: trend, seasonal and random. It’s easier to set out to find a suitable autoregressive and moving average (ARMA) model for the random part.
ts_1 <- ts(select2,frequency = 6)
de_1 <- decompose(ts_1)
plot(de_1)
We seek to fit a stationary Gaussian ARMA(p,q) model with parameter vector \(\theta=({\phi}_{1:p},{\psi}_{1:q},\mu,\sigma^2)\) given by \[ {\phi}(B)(Y_n-\mu) = {\psi}(B) \epsilon_n,\] where \[\begin{eqnarray} \mu &=& {\mathbb{E}}[Y_n] \\ {\phi}(x)&=&1-{\phi}_1 x-\dots -{\phi}_px^p, \\ {\psi}(x)&=&1+{\psi}_1 x+\dots +{\psi}_qx^q, \\ \epsilon_n&\sim&\mathrm{ iid }\, N[0,\sigma^2]. \end{eqnarray}\]
Let’s tabulate some AIC values for a range of different choices of \(p\) and \(q\). The ARMA(2,2) model provides a good preditve power. We will choose it for the subsequent complicated model.
MA0 MA1 MA2 MA3 MA4 MA5 AR0 2226.57 2209.52 2034.27 2032.62 1915.72 1905.42 AR1 2222.87 2205.10 2035.98 2028.27 1913.81 1902.18 AR2 2124.90 1899.59 1836.89 1838.53 1839.91 1841.32 AR3 2039.88 1858.61 1838.44 1840.57 1842.43 1843.29 AR4 1987.00 1848.36 1839.58 1842.37 1837.83 1844.04 Let’s fit the ARMA(2,2) model recommended by consideration of AIC.
arima_1<-arima(rand_1, order = c(2,0,2))
arima_1##
## Call:
## arima(x = rand_1, order = c(2, 0, 2))
##
## Coefficients:
## ar1 ar2 ma1 ma2 intercept
## 0.8745 -0.7877 -1.6479 0.6479 -9e-04
## s.e. 0.0396 0.0346 0.0576 0.0572 8e-04
##
## sigma^2 estimated as 13.2: log likelihood = -912.45, aic = 1836.89- In the first plot, residuals generally distributed aroud the line. Thus the residuals follow normal distribution. The Shapiro-Wilk normality test also shows we can not reject the null hypothesis. In the third plot, the residuals generally distributed aroud the line=0. In the fourth plot, residuals are generally inside the dashed line showing pointwise acceptance regions at the 95% confidence level under a null hypothesis of no correlation between noise. Therefore, we can not reject the null hypothesis of Gaussian noise.
##
## Shapiro-Wilk normality test
##
## data: arima_1$residuals
## W = 0.99509, p-value = 0.3684
5. SARMA Model
Let’s continue by adding the seasonable component. The SARMA\((2,2)\times(P,Q)_{6}\) model is \({\quad\quad\quad}{\phi}(B){\Phi}(B^{6}) (Y_n-\mu) = {\psi}(B){\Psi}(B^{6}) \epsilon_n\), where \(\{\epsilon_n\}\) is a white noise process and \[\begin{eqnarray} \mu &=& {\mathbb{E}}[Y_n] \\ {\phi}(x)&=&1-{\phi}_1 x-{\phi}_2x^2, \\ {\psi}(x)&=&1+{\psi}_1 x+{\psi}_2x^2, \\ {\Phi}(x)&=&1-{\Phi}_1 x-\dots -{\Phi}_px^P, \\ {\Psi}(x)&=&1+{\Psi}_1 x+\dots +{\Psi}_qx^Q. \end{eqnarray}\]
We need to tabulate some AIC values for a range of different choices of \(P\) and \(Q\). The \(SARMA(2,2)\times(2,1)_{6}\) model is chosen based on the AIC table.
SMA0 SMA1 SMA2 SMA3 SAR0 2231.19 2229.76 2227.34 2229.33 SAR1 2228.80 2227.12 2228.86 2230.44 SAR2 2226.62 2225.54 2227.07 2228.15 SAR3 2228.53 2227.47 2228.48 2227.65
sarima_1<-arima(select2,order = c(2,0,2),seasonal=list(order=c(2,0,0),period=6))
sarima_1##
## Call:
## arima(x = select2, order = c(2, 0, 2), seasonal = list(order = c(2, 0, 0), period = 6))
##
## Coefficients:
## ar1 ar2 ma1 ma2 sar1 sar2 intercept
## 0.7825 -0.8281 -0.5906 0.5632 0.1422 0.1296 1.5357
## s.e. 0.0796 0.0759 0.1103 0.1015 0.0662 0.0628 0.4305
##
## sigma^2 estimated as 38.87: log likelihood = -1105.31, aic = 2226.62
Therefore, the SARIMA model is acceptable.
6. Forecast
- To better examine the model constructed, 90% of the data is taken as train data and 10% as test data to do the forecasting. Most of the points were captured in the 95% confidence interval except some extreme peaks.
train <- select2[1:((0.9)*length(select2))]
test <- select2[(0.9*length(select2)+1):length(select2)]
train22 <- arima(train,order = c(2,0,2),seasonal=list(order=c(2,0,0),period=6))
pred_ <- predict(train22, n.ahead = (length(select2)-(0.9*length(select2))))
ts.plot(test,pred_$pred,col=1:2)
u=pred_$pred+1.96*pred_$se
l=pred_$pred-1.96*pred_$se
xx = c(time(u), rev(time(u)))
yy = c(l, rev(u))
polygon(xx, yy, border = 8, col = gray(.6, alpha = .2))
lines(pred_$pred, type="p", col=2)
7. Conclusion and Discussion
The \(SARMA(2,2)\times(2,0)_{6}\) model has been successfully constructed to simulate the MMC phase 3 motility pattern, though some further optimation might be needed to account for the extreme peaks;
The unsatisfactory part can be explained as the data has an unstable standard deviation and is right skewed. However, many time series analysis methods are built on the assumption of weak stationary: i. constant mean value function; and ii. the autocovariance function only depends on the length of time interval function;
- In the next step, we will
- aggrregate all the phase 3 manometry data of the subject in mid-duodenum throughout the 7hrs to get a better idea of the motility pattern along the time course
- accumulate the phase 3 manometry data of all the subjects at different sites along the gastrointestinal tract, build up a hierachical model to investigate the intra- and inter- subject variance and find out some significant covariable.
8. Reference
- DiMagno EP. Regulation of interdigestive gastrointestinal motility and secretion. Digestion 1997; 58 Suppl 1: 53-5.
- Schuster MM, Crowell MD, Koch KL. Schuster Atlas of Gastrointestinal Motility in Health and Disease. 2 ed. London: BC Decker Inc; 2002.
- Barrett KE. Gastrointestinal physiology. McGraw-Hill’s AccessMedicine. 2nd ed. New York, N.Y.: McGraw-Hill Education LLC,; 2014. p. xii, 321 p.
- Orlaith Burke, “Classical Decomposition” Statistical Methods Autocorrelation: Decomposition and Smoothing(2011): 15-16.