Let's imagine that the classifier plasmids are ready to be used in vivo, thus, in clinical trials.
- How can we predict the path, which the genetically engineered bacteria will follow inside the large intestine?
- How can we be certain that our bacteria will eventually reach the area around the cancerous tissue?
- Will our bacteria attach to the carcinoma or will the Colon Chyme wash them out?
- How do peristaltic movements of the colon affect the environment around the tumor (velocity field, pressure, etc)?
The purpose of our Colon Fluid Dynamics model is to give an answer to those questions and to further explore the unknown and chaotic world of the human large intestine.
In our project, the RNAi-based apoptotic circuit is to be eventually distributed via the gastrointestinal system towards the tumor inside the large intestine. Thus, it is essential that we examined carefully the properties of the colon's environment.
The large intestine is the last section of the gastrointestinal tract that is responsible for absorbing water and nutrients and concurrently converting digested food into feces. (Even though the length of the colon is considerably smaller than the small intestine's, the thickness of the walls as well as its diameter are bigger.) The length of the colon approaches the mean value of 1.5m and 6-7cm in diameter in most humans.
The Chyme (intestinal fluid containing the digested food) is transported from the small intestine towards the large intestine via the ileocecal sphincter. Chyme “settles” in the cecum, after its entrance into the large bowel, where it is mixed with bacteria which naturally reside in it constituting the gut's flora and contribute to much of the large intestine's functionality. The chyme is transferred consequently from one haustra to another by slow peristaltic waves. Those waves can be either mass movements, which are very slow and widespread movements happening a few times throughout the day and associated with food consumption by means of the gastrocolic and duodenocolic reflex or segmentation movements that mostly serve to chop and mix the intestinal content.
Peristalsis is a radially symmetrical contraction and relaxation of muscles that propagates as a wave down the tube (in our case the intestinal tube), in an anterograde direction. In the case of the human intestinal tract, smooth muscle tissue contracts in sequence to produce a peristaltic wave, which propels fecal aggregates along the tract. This is the main kind of movement that will be taken into account in our model. More technical details on the physical model of peristalsis are presented in the Human Colon Fluid Dynamics section of our model.
The large intestine can be considered a mechanical propulsion system implying that one should determine its elastic and viscous properties before proceeding with the model.
Some primitive tensile properties of the human large intestine are the maximal stress and destructive strain, valued at 0.9 MPa and 180% respectively[1]. In addition, to gauge the elastic properties of the intestinal wall one needs to determine the elastic modulus (Young's modulus) as well as the Poisson Ratio.
Young's modulus (E) is a numerical constant that describes the elastic properties of an elastic material undergoing tension or compression unidirectionally and functions as ameasure of the material's capacity to withstand alterations in its shape. The mathematical formula describing this physical entity is the following:
\[E = \frac{{stress}}{{strain}} = \frac{{F{L_0}}}{{A({L_n} - {L_0})}}\]The physical units of E in the metric system are Newtons per square meter (N/m2).
Poisson's Ratio (σ) is the ratio of transverse contraction strain to longitudinal extension strain in the direction of the stretching force applied on the material's surface. The large intestinal tissue is commonly considered as soft tissue and the mathematical formula describing its Poisson's ratio is the following:
\[\sigma = - \frac{{d{\varepsilon _{trans}}}}{{d{\varepsilon _{axial}}}}\]Those values exist in the literature; as a mean value across all locations of the colon (distal, medial and proximal) as well as in every possible orientation (circumferental and longtitudinal) for the Elastic Modulus equal to 5.18 MPa and as a commonly accepted value of 0.5 for incompressible materials for soft tissue, such as the colon.[2,3]
Colon cancer constitutes one of the most commonly diagnosed cancers worldwide, thus one could easily consider the importance of determining its morphology. More often than not, colorectal cancer emerges as a malignant transition of a polyp, therefore one could assume that they share their macroscopic geometry, particularly in the early stages.
In our Colon Cancer CFD model, we assumed that the polyp, which evolved to cancer, is found in the proximal ascending colon of the large intestine as it is one of the most common colon regions where carcinomas are found (together with the Sigmoid Colon). Thus, the geometry of our model represents the caecum and a major part of the ascending colon. The CAD file for the geometry was designed with the “Autodesk – Inventor Professional” software.
It contains 2D identical structures which were then revolved into 3D objects and unified as an assembly so that they form a 3D solid structure. The round surface on the left of the solid represents the fluid inlet, thus it could be considered as the final part of the caecum |
region close to the ileocecal sphincter. After the inlet, one could see the characteristic “curved cylindrical” shapes which represent the intestinal haustra. On the upper right side of the image is the fluid outlet which could be considered as the section close to the hepatic flexure which then continues to the transverse colon.
After the construction of the outer intestinal geometry, the geometry of the cancerous polyp was added to the structure. By adding a 2D sketch on a working plane parallel to one of the haustra and by constructing a mushroom-like shaped solid, we were able to design the tumor inside the colon.
The tumor is the blue object shown on the diagram. The solid which we created is very similar to the anatomical morphology of a human polyp. |
In this section we examine the physical laws behind the flow of the gastrointestinal fluids.
The basic differential equations that govern the flow of fluids are the Navier – Stokes equations and the continuity equation. The mathematical form of those equations is as follows:
\[\rho \left( {\frac{{\partial {\bf{V}}}}{{\partial t}} + {\bf{V}} \cdot \nabla {\bf{V}}} \right) = - \nabla p + \rho {\bf{g}} + \nabla \cdot \left( {\mu (\nabla {\bf{V}} + {{(\nabla {\bf{V}})}^T})} \right) - \frac{2}{3}\mu (\nabla \cdot {\bf{V}}){\bf{I}})\] \[\rho \nabla \cdot {\bf{V}} = 0\]where ρ (rho) stands for fluid density, V is the velocity vector of the fluid which produces a velocity field and μ (mu) stands for dynamic viscosity of the fluid. Let's analyze further each term in the equations:
- \[\rho \left( {\frac{{\partial {\bf{V}}}}{{\partial t}} + {\bf{V}} \cdot \nabla {\bf{V}}} \right)\] : This term represents all the inertial forces.
- \[ - \nabla p\] : This term all the pressure forces.
- \[\nabla \cdot \left( {\mu (\nabla {\bf{V}} + {{(\nabla {\bf{V}})}^T})} \right) - \frac{2}{3}\mu (\nabla \cdot {\bf{V}}){\bf{I}})\] : This term all the viscous forces.
- \[\rho {\bf{g}}\] : This term represents the external gravitational force of weight.
All these terms combined together, comprise the well-known Navier-Stokes which represent the conservation of momentum. The continuity equation for incompressible fluids represents the conservation of mass in the system. Another way of writing the Navier – Stokes equations for incompressible fluids is the following one:
\[\rho \left( {\frac{{\partial {\bf{V}}}}{{\partial t}} + {\bf{V}} \cdot \nabla {\bf{V}}} \right) = - \nabla p + \rho {\bf{g}} + \mu {\nabla ^2}{\bf{V}}\]These equations generally hold for incompressible Newtonian fluids. In the case of fecal matter, we consider the fluid to be Newtonian as it is composed mainly of liquid water (close to 90%) while inside the caecum and ascending colon and because the relationship between the viscous stresses and the local strain rate is linear.
Once we know the velocity field V as well as the pressure p in every spatial point of our geometry, we can determine the stress tensor, which in mathematical terms stands as follows:
\[{\bf{\sigma }} = - p{\bf{I}} + \eta [\nabla {\bf{V}} + {(\nabla {\bf{V}})^T}]\]From the stress tensor one could determine the values of two other important fluid mechanical quantities, the hydrodynamic force and the hydrodynamic torque, which are expressed mathematically as follows:
- Hydrodynamic Force: \[F(t) = \iint\limits_S(\[{{\bf{\sigma }} \cdot {\bf{n}}}\]) dS]
- Hydrodynamic Torque: \[L(t) = \iint\limits_S(\[{{\bf{x}} \times ({\bf{\sigma }} \cdot {\bf{n}})}\]) dS]
These two integrals are calculated on the surface S and x stands for the positions on that surface.
We should refer to a common dimensionless group in fluid dynamics that we used in our model. The ratio of inertial forces to viscous forces, Reynold's number (Re), is one of the most important dimensionless groups in fluid dynamics. The mathematical representation of Re is the following one:
\[{\mathop{\rm Re}\nolimits} = \frac{{\rho Vl}}{\mu }\]where ρ is the fluid density, V the fluid velocity, l the characteristic length (characteristic area in 3D objects) and μ the dynamic viscosity. The Reynold's number could be used as an indicator in order to determine whether a specific kind of fluid motion is laminar or turbulent of even transitional. In our case, the Reynold's number has a very low value equal to 0.00003 and thus, the viscous effects dominate the fluid. The reason why Reynold's number takes such a low value is that if we zoomed in the fluid we could see bacteria swimming at the same velocity as the fluid velocity field. Assuming that the size of a bacterium is about 10-6m, their swimming velocity is 30x10-6 m/s and the density of the fluid in which the bacteria swim is 1000 kg/m3 and the dynamic viscosity is around 0.001 Pa.s, the Reynold's number turns out to be equal to 1x10-5.
Life at low Reynold's numbers: In environments such as a fluid with a very low Reynold's number, the particles which “live” inside the fluid are considered “swimmers” as they deform their surface to sustain their movement. In our case, the E. coli bacteria, which constitute an important part of the gut's flora, have helical flagella to help them with their locomotion.
If a bacterium suddenly stops deforming its body, it will become a “victim” of the inertial forces of the fluid and thus, it will decelerate. At low Reynold's numbers, the drag force acting on the bacterium takes the form of viscous forces, \[{f_{drag}}\mathop \eta \limits^\~ UL\] , and the bacterium can reach a distance of about 0.1 nm before stopping its movement. In contrast, at high Reynold numbers the distance is longer and thus, Re at low values can be interpreted as a non-dimensional coasting distance.
One could describe the locomotion of a bacterium in low Re as a function of the swimming gait \[{u_s}(t)\] (velocity field on the body surface). At every instant, it can be assumed that the body is solid with unknown velocity \[U(t)\] and rotation rate \[\Omega (t)\]. Thus, the instantaneous velocity on the body surface is \[u = {\bf{U}} + {\bf{\Omega }} \times {\bf{x}} + {{\bf{u}}_s}\]. To calculate in every time step the values of \[U(t)\] and \[\Omega (t)\], one needs to determine both the velocity and stress fields in the problem in order to utilize the following integral:
\[{\bf{\hat F}} \cdot U + {\bf{\hat L}} \cdot \Omega = - \iint\limits_S(\[{{u_s} \cdot \hat \sigma \cdot n}\]) dS]Assuming that the swimmer is a bacterium with filament, what makes its movement in low Re possible is mainly the existence of drag forces perpendicular to the motion of the filament. The propulsive force generated along the filament (length L and deformation amplitude y(x,t)) is given by:
\[{{\bf{F}}_{prop}} \approx ({\xi _ \bot } - {\xi _\parallel })\int_0^L {\left( {\frac{{\partial y}}{{\partial t}}\frac{{\partial y}}{{\partial x}}} \right)dx{e_x}} \]This equation is a consequence of Purcell's Scallop Theorem and is derived for Re -> 0[4].
The software for the model's numerical simulations is Comsol Multiphysics (v. 5.0) thus, the numerical method utilized is Finite Element Method (FEM).
The number of physical problems that can be solved analytically is very small compared to all the real world physical systems that can be solved only by approximations. One way to find a solution to a complex physical problem such as fluid flow in the human large intestine is to utilize discretization methods that lead to numerical solutions. An example of space discretization method is the Finite Element Method which considers an analytical function as a linear combination of basis functions.
Comsol Multiphysics uses this method in order to compute the numerical solutions of fap in every single node of the space discretization. However, in order for such calculations to be possible one needs to discretize the problem geometry into nodes and elements. The grid that is being produced by the software user is called a “Mesh” and the shape of its elements can vary. In our model, the tetrahedron mesh was utilized, which is the 3D representation of the 2D triangle mesh. Although the meshing was implemented by the software itself by using “physics-based” mesh, the discretization method was tetrahedral as it is easy and quick to create but it also creates an unstructured grid which in our case is helpful as there are domain regions with large difference in the required scales. For instance, in the proximity of the tumor we needed a finer meshing as the numerical solutions would converge to more accurate values that way.
The mesh generated by Comsol Multiphysics is a physics-controlled mesh with finer elements. The tetrahedral space discretization is visible on the mesh plot. We present here two mesh plots; the first shows the tetrahedral meshing of the intestinal geometry and the second one depicts the interior region where the tumor meshing is visible. |
|
In CFD models one needs to determine the initial values of the parameters as well as the boundary conditions of the problem in order to get an accurate result from the solver.
The normal velocity of Chyme inside the small intestine is below 1 cm/s and inside the large intestine it's even slower. Thus, we assumed that a logical value for the inlet normal velocity is 0.5 cm/s. The outlet boundary condition is pressure and its value is 101.16 kPa, close to 1 atm [5].
The intestinal “tube” in our geometry is considered to be a no-slip wall in our study which implies that the velocity of the fluid particles close to the wall will tend to zero. The tumor region is a no-slip interior wall.
The velocity field of the fluid flow inside the large intestine is quite complex, when peristalsis occurs, as the intestinal walls contract in a periodical manner in order to forward the colon contents towards the rectum where they are being secreted. However, that kind of study is beyond the scope of our model. In our model the inlet velocity is a consequence of the contraction of the caecum but we assume that the time scale of the peristalsis is bigger than the fluid movement's; thus, no peristaltic waves occur throughout our model and the fluid is driven by its inlet velocity and a pressure outlet which is close to 1 atm.
Also, viscosity is a rather important factor when it comes to fluid flow. The intestinal fluid, due to nutrients and water absorption, doesn't possess a constant viscosity value; although the site of absorption and phase change of Chyme is the Transverse Colon and not the Caecum and Ascending Colon, the digested food determines to a high degree the success of nutrients and water absorption. Thus, we considered essential to examine the fluid flow under different dynamic viscosity values.
The dynamic viscosity in this section is the value of Water Dynamic Viscosity at T = 37oC. The following diagram shows the velocity field distribution inside the colon geometry.
The following diagram shows the velocity field distribution inside the colon geometry.
\[\mathop {\lim }\limits_{x \to \infty } x\]The dynamic viscosity in this section is the value of Water Dynamic Viscosity at T = 37oC. The following diagram shows the velocity field distribution inside the colon geometry.
In these diagrams, the magnitude of the fluid's velocity field is given with respect to the arc length of the tumor. In the following diagram, we depict the boundary arc chosen.
One could observe that between the velocity magnitude maxima appear a few local minima which we decided to call “Velocity Field Valley Traps (VFVT)” (or Valleys of the Muses). In these regions, the velocity around them has a greater value, thus, the bacteria “swimmers” that will reach these regions will most likely stay there and attach to the tumor rather than being washed out by the fluid. We could consider those regions as local stagnation regions of the fluid. Although the local maxima for each dynamic viscosity value are found at about the same location on the arc boundary of the tumor, the local minima or VFVTs don't follow the same rule. Having said that, the VFVTs for viscosity values 0.01 Pa.s and 1 Pa.s are located on the same regions which implies that for viscosity changes in that range the regions of the tumor on which the bacteria would attach doesn't change. Those regions are:
However, for viscosity value 0.6913 mPa.s the VFVTs seem to have a slight displacement from the other VFVT locations. More specifically, the first and third VFVTs “moved” towards the top arc length which implies that for low viscosities the VFVTs appear closer to more curved regions of the tumor geometry.
In our study of the fluid flow near the tumor, it is essential to determine the drag forces on the tumor if we want to predict whether our ‘bactofectors' will manage to stay attached to the tumor. Thus, due to the small cohesive forces between the bacteria and the tumor surface, it is safe enough to assume that the total stress forces on the tumor due to the fluid will apply to the bacteria as well.
We computed the total drag forces on the tumor arc boundary (same as the velocity magnitude) by integrating the total stresses over the entire arc length.
What seems interesting in these diagrams is the fact that for different values of fluid viscosity the total stress on the boundary remains the same. Why does this happen?
Well, let's remind ourselves that by integrating the total stresses on the arc boundary, we assume both the pressure and viscous forces. Thus, we should distinguish the two cases and examine which case is affected by the varying viscosity values.
If we could make a prediction, it would be logical to assume that viscosity affects the viscous forces which interest us.
