Applications of Computational Modelling and Simulation of Porous Medium in Tissue Engineering

In tissue engineering, porous biodegradable scaffolds are used as templates for regenerating required tissues. With the advances in computational tools, many modeling approaches have been considered. For example, fluid flow through porous medium can be modeled using the Brinkman equation where permeability of the porous medium has to be defined. In this review, we summarize various models recently reported for defining permeability and non-invasive pressure drop monitoring as a tool to validate dynamic changes in permeability. We also summarize some models used for scaffold degradation and integrating mass transport in the simulation.


Introduction
Tissue engineering refers to the in vitro regeneration or growth of a tissue or an organ for the purpose of biomedical studies in areas such as toxicology, tissue/organ repair, and medical device development.The technology is based on using biodegradable porous scaffolds to guide and support the in-growth of cells during tissue regeneration either at the site of grafting or in vitro [1].The scaffold transiently degrades leaving only the necessary healthy tissue.Currently, the primary goal of tissue engineering is to address the massive discrepancy between organ/tissue transplant needs and availability.Another goal is to use engineered tissue to develop safe and effective medication, treatments, and medical devices without threat to life [2].In vitro cell culture can help achieve these goals, although the growth of cells outside of the body can be difficult due to cell sensitivity.For this reason, cell culture environments closely mimic body conditions such as temperature, pH, and shear stress, as well as mechanical stimuli in some cases.Advancements in tissue engineering have led to the development of numerous techniques and housing environments in cell culture including bioreactors that aim to meet such requirements [3].However, due to the complex nature of cell systems, gaps exist in the knowledge of the interactions occurring within.
Tissue regeneration is a dynamic process where cell numbers are expected to increase.These cells also need to secrete matrix elements native to that tissue and matrix components need to assemble in the scaffold while the scaffold material degrades.These processes reduce the pore size available for fluid flow to reach a permeability matching that of the healthy tissues.Purely experimentally driven studies to determine parameter change and stimuli effects are costly in time and money, and do not provide a complete understanding of dynamic process [4].Computational modelling and simulation for optimization of in vitro tissue growth environments decreases experimental costs and allows easy variable manipulation to determine effect.Modeling is currently limited to activities that are capable of predicting or characterizing all of the processes that occur within a given cell culture system, but the success lies in the combination of well validated models into one simulation to tell the entire story of molecular metabolism, mass transport, scaffold decomposition, and many other phenomena on multiple levels.The combination of multiple models in simulation is necessary to build a complete understanding of what occurs throughout the tissue engineering process and what is needed to produce healthy, fully functional tissue [5].
The field of tissue engineering will find great success due to understanding gained from modeling and simulation compilations.This review will focus on the mathematical models currently available for describing mass transport in static and dynamic cell culture environments, as well as discuss methods for describing fluid flow effects on scaffolds.

Cell Culture
Numerous in vitro cell culturing techniques, such as growth on tissue culture plastic, cell suspension, cell aggregates, and bioreactors, have been investigated in an attempt to mimic body tissue.Each technique has advantages and faces varying challenges.Two-dimensional (2D) cell culture involves the growth of cells on pre-treated, flat, polystyrene plastic surfaces (Figure 1a).Cells attach to the stiff, plastic surface and spread until surface area runs out.Cells in 2D environments are submerged in nutrient-rich media, typically changed every 48 h, to promote cell vitality and proliferation.
Two-dimensional cell culturing techniques have had success in maintaining cell proliferation, and are easy to set up and maintain; however, 2D cultures are poor approximations of tissue in vivo.Two-dimensional cell cultures lack necessary architectural structure and face significant exposure limitations [6].Three-dimensional (3D) cell culture can provide necessary architectural structure, generally via a biocompatible, biodegradable scaffold.The scaffold is a porous medium that supports cell organization and motility and allows nutrient transport inside the 3D structure, creating an environment that more closely mimics tissue formations in vivo [7][8][9].In addition, 3D culture allows for higher cell density seeding over a given surface area, allowing for increased sensitivity and stronger signal strength during system analysis.Three-dimensional cultures are generally preferred to 2D cultures in tissue engineering, and for that purpose, the remainder of this review will focus on simulation of 3D cultures only.

Cell Culture under Static Conditions
Static cell culture systems are easy to maintain, but require regular media replenishment to maintain adequate nutrient supply.Mass transport of molecules, such as nutrients, can be modelled mathematically.The distribution of components in static condition is driven by diffusion as described by Fick's first law.Hence, theoretical models used in simulation apply some form of Fick's law to characterize mass transport in a system [10], where N i, di f f usive represents the diffusive mass flux through the system, D i represents the effective diffusivity of the component, ∇ represents the gradient operator used when evaluating a function in more than one dimension and ∇C i represents the gradient of concentration.The conservation of mass is given as, where t represents time.Substitution of Equation (1) into Equation (2) yields: Assumptions for simplification of a system are commonly made in simulation, especially in the beginning stages.For example, the assumption of a constant diffusion coefficient gives Equation (3) as: And for a steady state assumption, Equation (3) becomes: As most of the systems considered in tissue engineering evaluate consumption of nutrients or formation of signature products, a reactive term, R i , is added to Equation (3) giving [11]: The sign convention of the reaction term depends upon whether the molecule is being consumed (negative) or produced (positive).Unlike convection where nutrient transport can be estimated by the flow rate and concentration gradient, diffusion transport mainly requires evaluating the diffusivity of the nutrient molecule.The major factors affecting diffusion include porosity and morphology, scaffold biodegradation, swelling due to biodegradation, and shrinkage due to deposition of new biopolymers and cells.Standard approach to defining D i [12] in the literature of flow through porous medium is where τ is the tortuosity of the porous medium, ε p is the porosity of the scaffold, and D 8 is the free diffusivity of the component in water.Since determining τ of the porous medium is difficult, some use that of cells.Alternatively, Mackie-Meares relationship can be used to determine the effective diffusivity.An example of Mackie-Meares relationship [13] is Since there are many techniques to determine the porosity of the porous medium particularly in tissue engineering scaffolds, those porosity values can be directly utilized to determine effective diffusivity.

Cell Culture Involving Fluid Flow
Many 3D cell culture environments encounter challenges due to nonhomogeneous molecular distribution, particularly those on the interior of 3D cultures [14].Further, many body parts are exposed to stresses either due to the weight they carry (such as bone), function they perform (such as bladder and cartilage) or due to the flow of fluid (lung and blood vessels).These parts regenerate if mechanical stimulus is applied.In an effort to address distribution issues and applications of mechanical stimulus, some have investigated convection driven cell culture environments.In the case of bioreactors, convection can be introduced by a flow of fresh, nutrient-rich media into the environment as seen in parallel-flow and flow-through bioreactors.Convection can be introduced into a system via agitation or stirring, as seen in spinner flasks and rotating wall vessel bioreactors.Bioreactors of different flow configurations and sizes have been with an even wider array of nutrient distributing and mechanical simulating techniques [15].
Appropriate characterization of transport phenomena in cell culture simulation is vital to obtaining results that could be utilized for optimizing tissue regeneration.Mathematical modelling of fluid behavior, due to the addition of convection, increases the complexity and number of equations required to describe fluid behavior.Characterization of fluid behavior in a system where convection is present requires the addition of a momentum equation.Fluid momentum in non-porous regions is characterized by the Navier-Stokes equation written as: where terms respectively from Equation (10).In addition to these two assumptions, consideration of the fluid as incompressible with a constant density and viscosity would yield Equation ( 8) as: Use of the conservation of momentum equation must also include the use of the conservation of mass equation, or continuity equation, when solving.The continuity equation is given as: And if considered at steady state and constant density: Velocity selection does not affect equation selection for characterization of non-porous regions, but has a direct impact on the equation for characterizing the effect of fluid flow through porous scaffolds.In order to determine the applicable fluid flow equations, the Reynolds number within the pores (and/or Forchheimer number) is calculated.Based on the local Reynolds number, one can determine whether Darcy, Brinkman or Forchheimer equations are most appropriate.The Darcy equation is based on experimental observations of the one-dimensional flow of water at low velocities through a porous medium made of sand particles, and is written as: where X is the direction of fluid flow, v represents the superficial velocity, and κ represents the porous structure permeability.Many derivations for Darcy's law are used but one must be cautious, as this assumption is only valid when the local Reynolds number is below one [16].When the local Reynolds number is in the order of 100 or greater, the Forchheimer equation is used to incorporate the turbulent flow conditions.Forchheimer added a quadratic velocity term to the Darcy equation to account for microscopic inertial effects.The Forchheimer equation for one-dimensional flow is written as: where β is the non-Darcy coefficient.Taking into account local shearing effects between the fluid and porous walls, and the gradient pressure forces involved in transport of fluids from non-porous regions into porous regions, second order derivatives of the velocity are added to the Darcy equation, yielding the Brinkman equation: where X, Y, and Z represent directions that are mutually perpendicular.The Brinkman equation is useful when the local Reynolds numbers are in the range of one and 100.Brinkman equation can be viewed as an equation consisting of two terms (i) Darcy's equation and (ii) a viscous term, similar to Navier-Stokes equation, which provides molecular transport due to convection.
where N i, convective represents the convective mass flux, Ñ u represents the average velocity of the fluid, and again, C i represents the species concentration.Combining Equations ( 1) and ( 16) gives an expression for the transport of solutes within the fluid to give the advective-diffusive equation: Cell metabolism equations are not affected by the inclusion of convection for molecular transport.However, local concentrations of molecules are dictated by the advective-diffusive equation.Hence the rate of reactions are altered, which further leads to altered concentration profiles.

Incorporating Permeability
Properties of the porous scaffolds used in tissue engineering, such as biocompatible material, porosity, and tensile strength, play an important role in characterizing mass transport.The Darcy, Forchheimer, and Brinkman equations, which characterize fluid flow through porous scaffolds, require an important porous structure parameter, the permeability, κ.Sensitivity of model predictions are directly dependent on the accurate κ values, which depend on the geometric parameters of the pores in the scaffold [17].In order to predict flow properties with high fidelity, determining the permeability using physical characteristics of the porous scaffold such as porosity and the pore shape is an approach.Kozeny proposed a relationship in 1927 given as: where S represents the specific surface area, the ratio of total interstitial surface area to the bulk volume of the porous scaffold and K C represents the Kozeny constant, a dimensionless parameter that depends on the pore geometry.In human physiology, permeability of various tissues is defined using this Kozeny definition [18].The above expression of κ does not depend on fluid viscosity and density, which helps in extrapolating the simulation results to different fluids and flow conditions.There are a number of approaches available for calculating permeability, based on fiber orientation and size or pore area and number, depending on method and materials used for porous scaffold development.When scaffolds are formed through the process of electrospinning, fibers are randomly oriented (Figure 2a).
In order to calculate the permeability in randomly packed fibers, there have been many correlations.
One popular equation that seem to agree with the electrospun scaffolds is: κ " 3r 2 20 `1 ´ε p ˘`´ln `1 ´ε p ˘´0.931 ˘ (19) where r is the radius of the fibers.Since these scaffolds are very thin relative to their surface area, one could assume negligible porosity across the thickness.Then, porosity can be estimated using digital micrographs collected from scanning electron microscopy.Image analysis is performed to calculate the ratio of the open pore area to the total area of the image analyzed, which is considered as the porosity of the scaffold [19].If the fiber sizes are large, another correlation is developed using Lattice-Boltzmann method.In that case, the permeability is calculated using: Other methods do not use fibers to create scaffolds, but rather create pores within a solid medium.For example, freeze drying of water-based polymer solutions such as collagen, chitosan, and alginate produce pores mimicking water crystals.Permeability for such scaffolds are calculated [20], assuming pores to be cylindrical (Figure 2b), and calculated using the equation, where d is the pore diameter and n A represents the number of pores per unit area.However, some porous medium preparation techniques using leaching salt from organic solvent-based polymeric blocks such as polycaprolactone and poly lactic-co-glycolic acid (Figure 2c).One could assume the pores to be rectangular and then calculate the permeability using the expression: where L is the pore length, and W is the pore width.Utilizing permeability equations based on pore geometric characteristics such as pore size, shape, and n A is advantageous because experimental data can be used from scaffold analyses, typically reported in tissue engineering literature.Scaffold permeability depends on a combination of factors including pore size, porosity, pore geometry, and pore distribution [21].Alterations in any of these parameters affects permeability, and changes in scaffold permeability directly affects the structural integrity of the scaffold [22].Increased pore size decreases solid material volume and, as a result, the strength of the material.Scaffold strength is of particular importance during the initial cell seeding, as the extracellular matrix (ECM) that provides stability in vivo is not yet formed.In order to mimic the ECM and promote appropriate cellular behavior, scaffold strength should match in vivo tissue strength [23].Permeability also influences diffusion of nutrients into and waste out of the porous medium.The transport of nutrients necessary for cell survival into the scaffold, and waste metabolites which can be toxic to cells, is controlled by permeability [21].Thus, selection of appropriate scaffold permeability for achieving desired effects is critical.
The majority of scaffolds are designed with homogeneous porosity and pore distribution due to manufacturing and design limitations [24].Homogenous scaffolds or scaffolds of constant porosity are a common assumption employed in tissue engineering simulations.Many reports confirm that this to be a good assumption and useful in calculating scaffold permeability [25,26].However, many issues arise when considering uniform scaffold permeability.For example, natural tissue is not structurally homogeneous.Additionally, homogeneously porous scaffolds cannot address all of the mechanical and biological tissue requirements.Hence, heterogeneous scaffold development is currently of particular interest to researchers.Scaffolds with separate regions of porosity and scaffolds with porosity gradients have been investigated [27].In some cases, a multi-scale simulation approach is advantageous in accounting for environmental heterogeneity.Multi-scale approaches allow simultaneous investigation of the interaction and dependence of various factors on the micro (cellular) and macro (entire scaffold) scale [28].

Incorporating Scaffold Degradation
Characterization of scaffold degradation over time is important to producing and maintaining high-integrity, engineered tissue for in vivo use.If scaffolds degrade too quickly, the necessary structural architecture for promoting appropriate cell proliferation and functionality is lost, resulting in cell death.Computational modelling and simulation can provide a better understanding of the mechanism driving structural degradation.Some mathematical models consider scaffold degradation by chemical reactions and erosion due to fluid flow.Many studies have investigated the effects of environment (i.e., fluid velocity), porosity, and scaffold material on degradation rates by surface and bulk erosion, as well as by chemical reactions [29,30].For example [31], polymer degradation by hydrolysis is modeled using the following coupled equations, and where C e is the mole concentration of ester bonds of polymer chains, C m is the mole concentration of monomers, k 1 is the rate constant of the reaction which is independent of any catalyst, k 2 is the rate constant of a reaction catalyzed by the previously formed products, and γ is the disassociation power of the acid end group.These reactions are then related to the MW of the polymer present in the scaffold, which can be measured and then the model can be validated.The rate of scaffold degradation by hydrolysis is not restricted to the combination of Equations ( 23) and ( 24), but can also be characterized by an exponential decay or first order equation.Other research has considered enzymatic degradation according to Michaelis Menten kinetics in addition to degradation by hydrolysis [32].
Scaffold degradation can take on many forms and is highly dependent upon scaffold composition and component concentrations.

Incorporating Scaffold Deformation
Most porous scaffolds used in soft tissue regeneration i.e., except in bone regeneration, have low mechanical strength and are very elastic.In addition, porosity of the scaffolds affects their mechanical characteristics.Further, cells respond to stresses in the substrates [33].Hence, velocity selection in flow-through systems can directly impact cell vitality and scaffold structural stability.In order to understand how fluid flow deforms the porous scaffold, novel computational tools offer a moving mesh that can be used to predict deformation.The simulation needs generalized mechanical property of the scaffolds such as elastic modulus and Poisson ratio in drained conditions [20].Initial simulation with Navier-Stokes equation without moving mesh is used to create boundary condition such as the initial total load (F T ) applied by the fluid on the scaffold.This is calculated using, F T " ´n¨p´PI `τq (25) where n represents the normal vector.In addition to the stress due to the elastic nature of the scaffold, there is a hydrodynamic stress from fluid pore pressure, which is incorporated by coupling the strain using the relationship, where σ represents the Cauchy stress tensor, ε-represents the strain tensor, C represents the elasticity matrix, α B represents the Biot-Willis coefficient, and P f is the fluid pore pressure.Equation ( 26) reduces to Hooks Law of linear elasticity when α B is zero.Poroelasticity models in earth sciences for understanding soil consolidation use similar models in conjunction with Darcy's equation.However, these models have to use Brinkman or Forchheimer equations based on the flow regimes, as described above.Assuming the material to behave isotropically, and porosity to remain constant, elasticity matrix is obtained using the values of elastic modulus, E, and the shear modulus, G.In the linear region, shear modulus and elastic modulus can be coupled using the equation.
where ν represents the Poisson's ratio of the porous scaffold.Some have used quasistatic approximations and solved the equations analytically [34].Similar models have been attempted in cartilage and other tissues [35].In order to use these simulations, mechanical characteristics in particular elastic modulus and Poisson's ratio have to be analyzed in wet conditions similar to the condition where scaffolds are used.Currently, some of these properties are not available for many scaffolds and hence, they need to be determined.Nevertheless, adding these simulations help understand the local stresses experienced by cells.These can be then analyzed with the biological response of cells, and further optimization of the flow can be performed to improve the quality of the regenerated tissues.

Validation Techniques
One of the primary requirements for successful mathematical modelling and computational simulation is validation of the experimental results.Well-validated, mathematical characterizations of phenomena occurring within a tissue culture environment would help close knowledge gaps in tissue engineering.The computational approach summarized above utilizes scaffold properties routinely measured and reported in literature and pressure drop can be predicted by using an appropriate flow through porous medium equation.Further, pressure drop profiles for various flow regions can be simulated using a range of permeability encompassing scaffold permeability to healthy tissue.
In order to utilize these profiles, experimental validation can be performed at various flow rates using scaffolds alone.A simple flow-scheme for such purposes is shown in Figure 3. Placing a pressure transducer at the inlet of the bioreactor and another at the outlet will provide pressure drop across the scaffold.The pressure transducers would be connected to a computer to record pressure values for both the transducers.In addition, performing experiments without the scaffold could be used to account for pressure drop occurring without the scaffold.These experiments can help validate the simulations.The advantage of this approach is that pressure drop can be continuously monitored non-invasively, even when biological experiments are performed in sterile environment.The method allows for simulation validation without harming or contaminating the cell cultures.This is unlike many of the destructive assays performed for understanding the quality of the regenerated tissues.Permeability is expected to change during the course of tissue regeneration as cells will proliferate, synthesize extracellular matrix elements that are deposited in the porous scaffolds.Hence, the pressure drop increases as the permeability decreases.
In order to understand the implications of dynamic changes in permeability, simulations with decreasing pore sizes can be utilized.Matching with the simulated pressure drop to flow rate correlation can be utilized to understand the changes in permeability, quality of the regenerated tissue, and adjusting the flow rate in order to ensure nutrient distribution.There are a few commercialized products available to measure the permeability of the porous scaffolds by Darcy's Law that is used to monitor bone growth.The primary caveat is that pressure drop is a significant function of the flow pattern.Hence, for each bioreactor configuration, few profiles have to be generated.Nevertheless, these advances of monitoring pressure drop in real time would be invaluable to the tissue regeneration process.

Conclusion
Computational simulation of flow through porous medium is gaining significant attention in tissue engineering.As we validate developed correlations between pressure drop and scaffold pore architecture for wide variety of porous media used in tissue engineering, they can be utilized as a standard quality assessment method for regenerated tissues.Furthermore, simulation and modeling can guide preparation of better scaffolds and bioreactors for regenerating high quality tissues.These can be used in transplantation, drug screening, or understanding of disease progression.Modelling and simulation will allow ideas to be tested without waste of materials.The growth of an entire organ requires the inclusion of multiple cell types, as well as a vascular system.Addressing the physical and mechanical requirements of multiple cell types in a single environment is extremely difficult and too complex for current application.Perhaps, modeling will help develop individualize medicine.This is the future of tissue engineering.

Figure 1 .
Figure 1.Structural comparison of (a) Two-dimensional cell culture and (b) Three-dimensional cell culture.
represents fluid pressure, µ represents the fluid kinematic viscosity, ∇u represents the gradient of fluid velocity, ρ represents fluid mass density, I represents the identity matrix, and F represents the sum of the external forces (like gravity) applied to the fluid.The Navier-Stokes equation is used for Newtonian fluids, such as growth medium used in cell culture.Assumptions, such as steady state operation and negligible external forces, are often applied to simulation, eliminating the Bu Bt and F

Figure 2 .
Figure 2. Schematic of different pore configurations in porous scaffolds used in tissue engineering.These are in the flow direction.(a) Fibrous scaffolds formed by electrospinning; (b) Scaffolds with nearly circular pores when formed using freeze-drying; (c) Scaffolds with rectangular pores when formed using salt-leaching technique.

Figure 3 .
Figure 3. Flow-scheme to determine the pressure drop across the porous scaffold.