Construction of the Small Intestine Epithelial Cell Model Based on Molecular Dynamics Simulation and Preliminary Exploration of Drug Intestinal Absorption Prediction


 In this study, molecular dynamics simulation was applied to the construction of small intestinal epithelial cell membrane and prediction of drug absorption. First, we constructed a system of a small intestinal epithelial cell membrane that was close to the real proportion and investigated the effects of temperature, water layer thickness, and ionic strength on membrane properties to optimize environmental parameters. Next, three drugs with different absorptivity, including Ephedrine (EPH), Quercetin (QUE), and Baicalin (BAI), were selected as model drugs to study the ability of drugs through the membrane by the free diffusion and umbrella sampling simulation, and the drug permeation ability was characterized by the free diffusion coefficient D and free energy barrier (△G) in the processes. The results showed that the free diffusion coefficient D’ and △G’ orders of the three drugs were consistent with the classical experimental absorption order, indicating that these two parameters could be used to jointly characterize the membrane permeability of the drugs.


Introduction
Oral drug administration is the most common and convenient method for clinical administration, while the small intestine is the main place for oral drug absorption. Most drugs need to pass through the cell membrane of small intestinal cells and then reach the systemic blood circulation through passive diffusion or carrier-mediated active transport. Therefore, predicting the intestinal absorption of oral drugs is of great signi cance in evaluating the therapeutic effect of drugs [1].
The most effective way to evaluate the intestinal drug absorption is through the intestinal absorption model, which is currently commonly used includes cell model, animal model, and arti cial bio lm model.
Caco-2 cell model, the most commonly used cell model, is widely used because it is close to the absorption environment of the human body [2][3][4]. Considering that there are other transport carriers in Caco-2 cells, such as peptides, organic cation transporters, organic anion transporters, etc., it may lead to a lower level of drug absorption than the actual absorption level [5]. As for the animal model, in vivo and in situ animal models have the disadvantages of high cost, long period and complex processing of biological samples, and in vitro animal models have problems with poor accuracy of experimental data, which make experimental animal models not suitable for early drug intestinal absorption research, although they can best re ect the actual absorption of drugs [6][7]. With the advantages of high throughput, low cost, and high exibility, the parallel arti cial membrane in ltration model is an ideal in vitro model for preliminary drug absorption prediction through bio lm [8]. However, the phospholipids in this model are weakly bound together by intermolecular interactions, making them uid, deformed, misplaced, and di cult to manipulate [9]. Therefore, the traditional drug absorption model is not suitable for the high-throughput screening of early drug intestinal absorption studies, and an e cient and simple model is urgently needed to evaluate drug absorption.
Molecular dynamics (MD) simulations, a new technology developed in recent years, expounds the mechanism of drug and bio lm action from micro and mesoscopic scales to obtain important information that is di cult to obtain in the experiment, so it is widely used in drug absorption [10][11][12]. Mojumdar etc. [13] applied it to the study of transmembrane absorption of articaine, a local anesthetic, and concluded that the effect of drug molecules on membrane electrostatic potential might be one of the molecular mechanisms of anesthetic effect. Boggara etc. [14], Simulated the diffusion process of aspirin and ibuprofen in DPPC membrane and compared the in uence of drugs with different charged forms on membrane permeability under the in uence of pH. Thus, it was concluded that the electrically neutral drug molecules passed through the phospholipid membrane through passive diffusion, while the electrically charged drug molecules might pass through the transient pore effect. The results above show that MD simulations study the effect of drugs on bio lm structure and the process of drug absorption through the membrane by establishing a cell membrane model from the point of intermolecular interaction. Therefore, it provides a new method for predicting intestinal drug absorption and the study of the interaction law between drugs and bio lm.
In this paper, the small intestinal epithelial cell membrane model was built based on the CHARMM-GUI platform. The effects of temperature, water layer thickness, and ionic strength on the model were investigated through MD simulation so as to optimize the parameters of the established model and make it closer to the real small intestinal epithelial cell membrane. In addition, the free diffusion and umbrella sampling processes of drugs with different absorption rates were simulated, and the parameters in the simulated membrane permeation process were used as evaluation indexes. By comparing the simulation results with the experimental results, the reliability of the model was proved, and the parameter evaluation method of intestinal drug absorption matching with the model was established in order to realize the rapid and simple prediction of intestinal absorption.

Optimization of model parameters
The performance of the small intestinal epithelial cell membrane was affected by temperature [18][19], water layer thickness, and ionic strength [20][21], respectively. The properties of the mixed membrane were studied by observing the changes of the above three factors so as to optimize the environmental parameters and establish the small intestinal epithelial cell membrane model close to the real composition. In this way, we can better predict the intestinal absorption of the drug. In this paper, Area per lipid (APL), Bilayer thickness, and diffusion coe cient D were selected as evaluation parameters to screen out the optimal model parameters.
The annealing simulation was used to study the phase transformation process of the mixed membrane, to investigate the in uence of temperature on the membrane properties. The annealing simulation time was 1.6 μs, and the whole temperature span was 160 K (230-390 K). For a bilayer system composed of only lipid molecules, the water layer thickness was generally 22.5Å, and the ionic strength was 0.15mol.
At this time, the lipid membrane can maintain the equilibrium state of water and ions [22]. However, the bilayer system constructed in this paper was composed of multiple lipid molecules. In order to screen out the optimal environmental parameters to meet the morphological requirements and equilibrium state of the lipid membrane, the effects of different water layer thickness and ionic strength were selected to investigate the established model. The water layer thickness was 10.5, 22.5, and 33 45 Å, and the ionic strength was 0.05, 0.15, and 0.25 mol, respectively.

Predictive characterization of drug absorption through the membrane
In this paper, three drugs, EPH, QUE, and BAI, previously studied by our research group, were selected as model drugs. The absorptivity experimental data of each drug was shown in Table 1, where -LogP was obtained by converting the apparent permeability coe cient (Papp) in the Caco-2 model. The lower the -LogP of a drug, the better the drug was absorbed, so the cell experiment results showed that the order of absorption of the three drugs was FER QUE BAI.
In the optimized small intestinal epithelial cell membrane system, we simulated the free diffusion of three model drugs into the cell membrane and the process of umbrella sampling through the cell membrane and selected two parameters, diffusion coe cient D and △G, to characterize the free diffusion rate and membrane permeability of drugs. By comparing the simulation results with the cell experimental results, the reliability of the small intestinal epithelial cell model was veri ed, and the parameter evaluation method for the prediction of oral drug absorption through the membrane was established at the simulation level.
To simulate the free diffusion of drugs into the cell membrane, drug molecules with a certain molar concentration were rst placed in the water layer of the cell membrane to replace the corresponding number of water molecules (Fig. 2). In order to simulate the process of drug molecules through the membrane, a drug molecule was placed in the center of the double lipid membrane, and the double lipid membrane model and each drug molecule were located in the origin coordinate (X = 0, Y = 0, Z = 0), and then by combining their atomic coordinates to form a new multi-component complex. Using umbrella sampling to apply a traction force to the drug molecule makes it pass through the lipophilic center and hydrophilic center of the lipid bilayer along the Z-axis until it reached the water layer outside the membrane so as to obtain the complete trajectory of the drug molecule passing through the membrane [23].

Evaluation parameters
APL, the area of a single lipid molecule in the lipid bilayer cross-section, can be used to illustrate the density of the lateral arrangement of small intestinal epithelial cell membrane molecules in the system and evaluate whether the system is in equilibrium according to its change over time. APL can be calculated from gmx energy in GROMACS.
Bilayer thickness, the thickness of lipid bilayer of the cell membrane, can be used to characterize the longitudinal stacking of cell membrane molecules in the small intestine epithelium. gmx density in GROMACS can be used to obtain the density of the lipid head group, and the bilayer thickness is obtained by calculating the distance of the peak in the mass distribution map of the lipid head group along the Zaxis.
Mean square displacement(MSD), the square of the distance from the initial position of the particle at t time. According to the curve of the MSD with time, the particle diffusion coe cient D can be obtained. For the lipid tail chain, the transverse diffusion coe cient D can characterize the uidity of the tail chain, while for the drug molecule, the Z-axis diffusion coe cient D can characterize the rate of its entry into the lipid bilayer. The diffusion coe cient D of drugs and phospholipid molecules can be calculated by gmx msd in GROMACS.
△G, the free energy change of drugs passing through bio lms, can investigate the di culty of different drugs passing through bio lms, and the smaller the △G is, the easier the drug is to penetrate the membrane. gmx wham in GROMACS can be used to transform the probability distribution of drugs in lipid membrane to obtain the free potential energy surface, through which the △G in the membrane penetration process of each drug can be obtained.

The results of model parameter optimization
The small intestinal epithelial cell membrane model preliminarily established in this paper is shown in Fig.3. In order to better predict the intestinal absorption of drugs, the in uence of temperature, water layer thickness, and ionic strength on the small intestinal epithelial cell model was investigated in this paper to select the optimal environmental parameters.

Effects of temperature on membrane properties
The changes of APL of AP and BL with temperature are shown in Fig.4 and Fig. 5. It can be seen from the gure that in the temperature range of 230-390 K, the in uence of temperature on both sides of the membrane is roughly the same, and the APL increases with the increase of temperature. For the heating process, there are two sharp points ( Fig.4a and Fig.5a). A smaller mutation point is between 280-290K, before which APL is small and almost unchanged, indicating that the lipid head groups are well arranged and tightly packed, and the membrane is in a complete gel state. While APL increased signi cantly after this mutation point, indicating that the arrangement of lipid head groups was looser than before, which was speculated to be related to the pre-phase transition of the cell membrane. The other mutation point was between 310-330K. After this mutation point, APL increased rapidly and did not appear saturation point with the increase of temperature. This indicates that with the increase of temperature, the thermal energy destroys the interaction between the lipid head groups, which increases the free space of the lipid head groups, leading to the complete transition to the liquid state. For the cooling process, APL decreases uniformly with the decrease of temperature, and there is no obvious mutation point on both sides of the membrane with the decrease of temperature. ( Fig.4b and Fig.5b).
The membrane thickness of AP and BL varies with temperature, as shown in Fig. 6. In the temperature range of 230-390 K, the in uence of temperature on both sides of the membrane is roughly the same, and the membrane thickness decreases with the increase of temperature. Moreover, for the heating process, there are corresponding turning points around the two temperatures at which the phase transition occurs, which is consistent with the results of APL. For the cooling process, there is no obvious turning point except that the AP lm has a big uctuation at 250K, which is also roughly consistent with the results of APL without mutation during the cooling process. In addition, to further explain the in uence of temperature on the membrane morphology, AP membrane was taken as an example to analyze the lipid molecular arrangement at three temperature points where the thickness of the membrane was signi cantly different during the heating process (Fig. 7). The results show that the higher the temperature is, the wider the peak shape of the density distribution map of the lipid head-based beads along the Z-axis is, which indicates that the higher the temperature, the more disordered the lipid molecules are arranged. On the contrary, the lower the temperature, the narrower the peak shape of the density distribution map along the Z-axis of the lipid head-based beads, indicating the more orderly the arrangement and the larger the peak distance between the two peaks, indicating the larger the lm thickness. This is also consistent with the results of changes in APL. The increase of temperature would destroy the interaction force between the lipid head-based beads, making them move more freely and arrange more disorderly.
The morphology of AP and BL varies with temperature, as shown in Fig. 8. As can be seen from the gure, at 250K and 280K, the lipid membrane is tightly arranged and orderly, and the membrane thickness is basically unchanged, so the membrane is in a complete gel state at this time. At 310K, the arrangement of beads was obviously looser than before, indicating that a pre-phase transition occurred at some time point between 280-310K, and the lm had already transited from gel state to liquid state, but the lm still maintained a good lm shape without severe deformation. When the temperature reached 340K, the lipid bilayer began to bend and deform, and with the increase of temperature, the lm morphology became more curly until the deformation was severe and no longer existed in the form of the bilayer, indicating that a phase transition occurred at some time point between 310K-340K. The above changes in membrane morphology with temperature more intuitively prove the previous results of APL and membrane thickness.
The above results indicated that the pre-phase transformation of the small intestine epithelial cell model was between 280K and 290K, and the phase transformation was between 310K and 330K, which was consistent with the reported cell membrane phase transformation temperature near the physiological temperature [24][25][26]. Moreover, it is worth noting that when the temperature is maintained between 290 K and 310 K, the gelatinous state of the membrane coexists with the liquid state, and the lipid bilayer of the membrane is in good shape without serious bending deformation, which is close to the real state of the membrane.

Effects of water layer thickness and ionic strength on membrane properties
The results of APL and membrane thickness changing with water layer thickness and ionic strength are shown in Fig. 9. When the thickness of the water layer is 22.5Å, APL is 0.46 nm 2 , and membrane thickness is 4.1nm. With the thickening of the water layer, the APL increases, and the membrane thickness also becomes thicker. This indicates that increased hydration causes more water molecules to contact with the lipid head region, which reduces the compactness of lipid molecules horizontally, thus increasing APL, while the entry of excess water molecules also increases the membrane thickness when the ionic strength increases or decreases, the changes of APL and membrane thickness tend to be stable, indicating that the ionic strength has little effect on the membrane morphology in the model built in this paper.
In conclusion, when the temperature is maintained between 290 K and 310 K, the APL of the membrane constructed in this paper is stable at about 0.46-0.47 nm. At this time, the lipid bilayer of the membrane is in good shape, and the membrane is neither in absolute gel state nor will there be severe bending deformation, but close to the real membrane state, which will be conducive to the follow-up study of drug permeation membrane absorption. Therefore, the simulated temperature we selected was 310 K, which was also consistent with the temperature of the cell membrane under physiological state. The water layer thickness was 22.5Å, and the ionic strength was 0.15mol, which was consistent with the relevant literature reports.

3.2Predictive characterization of drug absorption through the membrane
In this section, three model drugs QUE, EPH, and BAI, were used to simulate the complete transmembrane process of the drug by free diffusion and umbrella sampling in the optimized model. Free diffusion simulation, in which the drugs enter the cell membrane through free diffusion, can be used to re ect the di culty of drug penetration preliminarily. Considering the interaction and absorption process between drugs and cell membranes in the physiological environment, the longitudinal drug diffusion along the Z-axis is closer to the real drug permeation process. While the transverse drug diffusion along the X-axis mainly indicates that the drug has different degrees of disturbance to the lateral distribution of lipid molecules. Therefore, the MSD and diffusion coe cient D of the drug along the Z-axis were calculated respectively to characterize the free diffusion velocity of the drug, and the results were shown in Fig.10 and 11. According to the gures, although the proportion of lipid molecules on both sides of the membrane was different, the longitudinal diffusion coe cient D of each model drug also had some differences, but the overall diffusion trend remained the same, and the diffusion coe cient D was ordered as EPH > QUE > BAI. The results showed that the changes of lipid membrane components in the selected lipid components in this study had little effect on the drug diffusion rate, and the order of three drugs permeating the membrane re ected by free diffusion simulation was EPH > QUE > BAI, which was consistent with the results of cell experiment.
Umbrella sampling simulation, in which the drug passes through the intact cell membrane under the action of the traction force, can be used to re ect the drug penetration capacity. In the simulation, the free energy potential surface (Fig. 12) of each model drug through the membrane was obtained, and △G of each model drug through the membrane was calculated by using PMF graph data ( Table 2).
It can be seen from the table that the △G of each model drug through the membrane of both sides is different to some extent, but the same trend on both membranes of the △G re ected the same tendency of membrane permeability of each model drug, indicating that the change of lipid membrane composition in the selected lipid components in this study had little effect on the membrane permeability of the drug. In addition, it can also be seen from the table that the △G of each model drug is greatly different, and the descending order is EPH, QUE, and BAI, indicating that the membrane permeability of these three model drugs is EPH > QUE> BAI, which is also consistent with the experimental absorption order.
In this part, the free diffusion simulation, which preliminarily characterized the di culty of drug membrane penetration, combined with the umbrella sampling simulation, which evaluated the ability of drug membrane penetration, was used to re ect the ability of drug absorption. The results showed that the free diffusion coe cient D and △G order of the three drugs in two simulations were consistent, and the order was EPH > QUE > BAI, which is also consistent with the order of absorption in cell experiments. This indicates that the model established in this paper is reliable, and the two parameters can be used to re ect the absorption rate of drugs to a certain extent jointly.

Conclusion
In this study, we have built a model of the small intestinal epithelial cell membrane that is close to the real proportion and investigated the in uence of temperature, water layer thickness, and other factors on membrane properties to optimize environmental parameters. The results show that when the simulated temperature is 310 K, the water layer thickness is 22.5 Å, and the ionic strength is 0.15 mol, the membrane is in good condition, which is close to the real cell state.
We also simulated the process of free diffusion and umbrella sampling of different drug models and characterized the drug permeation ability by the changes of parameters during the process of drug molecules permeating the membrane. The results show that the free diffusion coe cient D and free energy barrier of the model drugs re ect the same tendency as the classical experimental results, which can jointly characterize the absorption rate of the drug.
The focus of our work is to build a model of small intestinal epithelial cells close to the human body and to establish a parameter evaluation method for the prediction of oral drug permeable membrane absorption at the simulation level. Compared with the traditional intestinal absorption model, this model has the advantages of saving time, resources and manpower, and is suitable for high throughput screening for early studies of intestinal drug absorption. Future work will focus on the study of the drug permeation mechanism Declarations Funding This work was funded by the National Natural Science Foundation of China (Grant No.81473364).

Con icts of interest/Competing interests
The authors declare that there is no con ict of interests regarding the publication of this article.

Availability of data and material
The datasets used or analyzed during the current study are available from the corresponding author on reasonable request.

Code availability
The experiment was carried out on the basis of Gromacs 5.0.7, Packmol 18.013, Visual Molecular Dynamics 1.9.2, and CHARMM-GUI.
Authors' contributions YS and XS contributed to the conception of the study; YS, MS and XD performed the experiment; QQ, YL, and LL participated in the data collection and the analysis of the results; YS and XS aided in drafting and revising the manuscript; JY, ZL, and QZ helped perform the analysis with constructive discussions.

Consent to participate
Using the small intestinal epithelial cell membrane model established in this study and three representative drugs as models, the simulated force eld was fully validated.

Consent for publication
This study provides a model of the small intestinal epithelial cells close to the true composition and a new method for predicting intestinal drug absorption.   Density distribution diagram of lipid head beads along the Z-axis during the heating process of AP membrane. Black is 230K-250K, red is 310K-330K, and blue is 370K-390K. Effects of water layer thickness (left) and ionic strength (right) on membrane properties. Red represents lm thickness, and blue represents APL. MSD of each drug in the z-axis direction of AP and BP membranes. A: AP-Z; B: BL-z. Red is BAI, blue is EPH, and green is QUE.
Page 18/19 Figure 11 diffusion coe cient D of each drug in the z-axis direction of AP and BP membranes. Red is BAI, blue is EPH, and green is QUE. Figure 12