Immersed Boundary Models for Quantifying Flow-Induced Mechanical Stimuli on Stem Cells Seeded on 3D Scaffolds in Perfusion Bioreactors

Perfusion bioreactors regulate flow conditions in order to provide cells with oxygen, nutrients and flow-associated mechanical stimuli. Locally, these flow conditions can vary depending on the scaffold geometry, cellular confluency and amount of extra cellular matrix deposition. In this study, a novel application of the immersed boundary method was introduced in order to represent a detailed deformable cell attached to a 3D scaffold inside a perfusion bioreactor and exposed to microscopic flow. The immersed boundary model permits the prediction of mechanical effects of the local flow conditions on the cell. Incorporating stiffness values measured with atomic force microscopy and micro-flow boundary conditions obtained from computational fluid dynamics simulations on the entire scaffold, we compared cell deformation, cortical tension, normal and shear pressure between different cell shapes and locations. We observed a large effect of the precise cell location on the local shear stress and we predicted flow-induced cortical tensions in the order of 5 pN/μm, at the lower end of the range reported in literature. The proposed method provides an interesting tool to study perfusion bioreactors processes down to the level of the individual cell’s micro-environment, which can further aid in the achievement of robust bioprocess control for regenerative medicine applications.


Introduction
The culture of stem cell populations in dynamic set-ups, for example in perfusion bioreactors, holds great potential for the production of tissue engineered constructs [1]. The use of bioreactors permits automated seeding and expansion of progenitor cells [2,3], facilitating the production of clinically relevant cell populations in close systems, while maintaining their phenotype and bone forming potential [4]. Additionally, these systems can provide controlled biomechanical stimuli, such as fluid flow-induced shear stresses, that might significantly affect stem cell properties during dynamic culture in bioreactors. For example, mechanical stimuli have been associated to early stem cell lineage commitment [5] and osteogenic priming in the absence of inductive growth factors [6][7][8]. Moreover, they have been shown to further promote osteogenic differentiation of bone marrow, periosteum and adipose derived osteochondroprogenitor cells in the presence of osteoinductive growth factors [9][10][11][12][13]. Osteogenic differentiation has been linked to the magnitude of shear stress, showing dose dependent enhancement of extra-cellular matrix deposition and subsequent mineralization by the cultured cells [14][15][16][17].
In order to characterize the dynamic environment throughout cell seeded scaffolds in perfusion bioreactors, many Computational Fluid Dynamics (CFD) modeling studies have been presented in the past decade [18][19][20][21][22]. However, the majority of these studies considered empty scaffold geometries without incorporating a cell domain. Recently this issue was addressed by representing the growing neotissue as a porous medium in order to model the effect of neo-tissue growth on the flow profile [23][24][25][26][27]. Still, these models predict the local distribution of shear stress and pressure throughout a volume averaged porous domain and do not take into account the local mechanical and geometrical environment of individual cells.
Mechano-transduction of stress induced by shear flow conditions is highly localized at specific areas of the cell's interface with its environment, such as focal adhesions, FAs [28], and primary cilia [5]. The latter have been shown to be involved in the osteogenic response of bone cells to dynamic shear flow conditions [29], as well as in remodeling of the extracellular matrix [30]. The amount of force perceived at the level of FAs as a result of external flow conditions is influenced by the cell's mechanical properties, cell shape and the geometry of its microscopic environment, e.g. location of attachment points, and presence of extracellular matrix (ECM). In this respect, the concept of cell cortical tension has gained a renewed interest in the last years as a mediator of mechano-transduction processes [31]. Cortical tension is created by the cell itself through active acto-myosin contractility, resulting in a prestressed cytoskeleton. This self-generated stress is an essential aspect of the tensegrity theory as introduced by [32], which posits that the cytoskeleton constitutes a tensegrity structure, with tension generating cortical stress fibers as 'ropes' and load bearing capacity provided by other cytoskeletal elements, the This work is part of Prometheus, the Leuven R&D division of Skeletal Tissue Engineering. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
substrate or the ECM. Next to additional structural integrity, the mechanically stressed state of a cell can boost the mechano-sensitivity of a cell [33]. External flow however, can also contribute to locally elevated levels of cortical tension, especially close to attachment point such as FAs. This passive source of cortical tension, as well as its importance relative to the cell-generated active tension and prestress, has not yet been investigated for perfusion cell culture systems. Therefore, a scaling gap exists between small scale i.e. 'single cell' and 'neotissue/whole scaffold' macro-scale that needs to be bridged. Computational models using realistic single cell geometries are a prime candidate for facilitating this task. Similar concepts of cross-scale model integration in order to capture mechanical interactions across scales from a whole organ to single cell level have been already described and envisaged [34,35] and have served as a paradigm for the current study focusing mostly on a scaffold-based in vitro process.
Computational models of cell deformation due to shear flow have been developed considering the cell as a 2D Gaussian interface [36] or a 3D linear elastic solid [23, [37][38][39][40][41][42][43][44][45][46][47]. The latter use a mixed Lagrangian-Eulerian formulation to solve the Fluid-Structure Interaction (FSI) problem, with a coupling through continuity boundary conditions. Additional numerical methods have been recently developed for modeling fluid-flow driven solid deformations in a biomechanical context. Immersed finite element methods have been used for modeling soft tissue deformation under the influence of blood flow [47] and within the walls of the aortic root [48]. In addition cell motility and deformation through contracted channels reminiscent of microfluidic experiments were also captured using a similar method operating with a single analysis mesh for solid and fluid that was not subjected to any deformation [49]. For larger deformations, the interaction between cell and fluid has been resolved by means of the levelset method [50]. Alternatively, the Immersed Boundary Method (IBM) is able to explicitly take into account discrete entities in the cell's cortex and, possibly, its internal cytoskeletal structure. It has been used to model the movement and deformation of vesicles, red blood cells and bacteria under flow conditions [51,52]. An FSI model for osteoblasts attached to scaffold struts was recently published [53], with a rigid single cell consisting of a half-sphere with two focal adhesion points. In the work presented in this study, more realistic cell shapes are introduced, which are not rigid but deform due to the fluid flow. Still, the cytoskeleton constitutes a highly complex, mechanoadaptive material [54][55][56] and its mechanical behavior differs between various temporal and spatial scales, [57,58]. Hence at present, only a strongly simplified mechanical representation of a complete attached cell is considered computationally feasible.
The main purpose of this study is to use the IBM to investigate fluid-induced mechanical stimuli on progenitor cells used for bone tissue engineering (human periosteal derived cells, hPDCs) attached to regular pore titanium scaffolds inside a perfusion bioreactor set-up. Each cell is represented by a simplified model of the cortical shell, similar to [59], supplemented with discrete Focal Adhesions (FAs) and an elastic nucleus. A multi-scale modeling approach is presented, consisting of a CFD analysis at the scaffold macroscopic (tissue) scale in order to determine suitable input boundary conditions at the microscopic scale (single cell scale) where the fluid-structure interaction is modeled by means of the IBM. The impact of the spatial location of the cells within the scaffold during flow perfusion on a number of key mechanical quantities at the cellular scale was investigated. To illustrate how (location-induced) geometrical differences might affect the biomechanical environment of single cells, three characteristic locations and corresponding cell geometries were chosen: one cell spread along the direction of the flow (A), one facing the flow (F) and one bridging between two struts (B). Furthermore, the presented model was used to assess how small clusters of cells attached to the scaffold are mechanically affected by perfusion flow. In order to investigate mechanical effects of mutual shielding [60], a 'three cells configuration' (T) facing the flow was investigated.

Immersed Boundary Method
The IBM has been developed for simulating moving, deformable membranes immersed in a fluid, based on a combination of an Eulerian and a Lagrangian approach [61]. The deformable object (the cell in this study), is represented by a discretized membrane/cortex Γ(t) and is able to move freely through the fixed Eulerian mesh O on which the flow is computed. The interconnection between both lattices is accomplished by means of a smoothed Dirac function δ.
In the 3D mesh O, the equations for incompressible Stokes flow are solved (as appropriate for the low Reynolds numbers typically encountered in bioreactors, see also in the supplementary information): with suitable boundary conditions which are explained in following section. In Eqs (1) and (2), u represents the fluid velocity, p the pressure and μ the viscosity. The influence of the cell boundary Γ(t) immersed in the fluid is taken into account through the distributed force density F and can be expressed as: Here, x are the Eulerian coordinates and X(s,t) are the discretized cell membrane coordinates indicating the position of the membrane at time t. As mentioned previously, the interaction between both meshes is realized through the introduction of a Dirac function δ defined by the following continuous function of the distance r: Using Eq (4), Eq (3) can be rewritten in a discrete formulation with N the number of nodes of the cell membrane, h the Eulerian mesh size and Once the flow is computed, the membrane positions are updated using the following equation of motion: Mechanical representation of a cell The mechanical representation of a cell in the Eulerian domain relies on the work presented in [62] where the underlying mechanisms and assumptions are discussed in detail. Briefly, the model assumes that most of the cytoskeletal material is present in a relatively thin cortical shell and that an elastic description of deformations at short timescales is adequate. Hence, the mechanical properties of this cortical shell will account for the mechanical response of the complete cytoskeleton. The immersed boundary which represents the cell is composed of a triangulated surface with a stretching stiffness k s and a bending energy k b . Finally, the cell nucleus is represented as a submerged solid elastic sphere with Young's modulus E n .-see Fig 1. The linear spring force of node i for each connected node j at distance d ij and resting length d ij 0 is expressed as: where e ij denotes a unit vector pointing from i to j. A moment of bending is computed between all adjacent triangles k and m with angle θ km and resting angle y km 0 : A force corresponding to this moment is applied to the non-common points of each of the two triangles, and a compensating force is applied to the common edge points, ensuring that the total force on the cell remains unchanged, i.e. for two triangles with common nodes c 1 and c 2 and non-common nodes l k and l m : doi:10.1371/journal.pcbi.1005108.g001 L k and L m are the distances from resp. node l k and l m to the line containing the common edge, and n k and n m are the normal unit vectors of triangle k and m. We denote the total bending force contribution of all adjacent triangle pairs of node i as f i b . This type of bending stiffness is commonly found in the literature for Red Blood Cell models [63]. The cell's volume is maintained through an effective bulk modulus K. For this, an internal pressure P v is computed based on a cell's volume V and equilibrium volume V 0 : Subsequently, a force f i v ¼ P v A i n i is obtained for each node i with A i and n i respectively the area and outward normal unit vector of each node, both of which are calculated using a discrete version of the Laplace-Beltrami operator. Furthermore, the nucleus is represented as a solid, elastic sphere with Young's modulus E n for which contact with the cortical nodes is considered Hertzian, i.e.
for node i indenting a nucleus with radius R n and overlap distance d i n , and with e i n a unit vector pointing from the nucleus center to node i. Discrete attachment points serve as Focal Adhesions (FAs). These points are placed outside of the fluid domain and are therefore not displaced by the fluid. Finally, the total force per node f i (s,t) is computed as the sum of all aforementioned partial forces: Atomic Force Microscopy measurements of cell cortical stiffness Cell cortical stiffness was measured using Atomic Force Microscopy (AFM). Measurements were performed using a Nanowizard 3 BioScience AFM (JPK) with a working range of 100×100×15 μm mounted on the stage of an inverted microscope (Olympus 1) placed on a vibration-isolation table. A V-shaped gold-coated silicon nitride cantilever with a four-sided pyramidal tip (Budget Sensors) with a nominal tip radius r tip of 15 nm and an opening angle θ of 35 degrees was used as the probe. The spring constant k spring of the cantilever was ca. 0.3 Nm −1 . Exact values have been calibrated using the thermal fluctuation method. Force curves have been recorded at 5 μm/s approach and retract speed, of which only the approach curves have been analyzed to arrive at the instantaneous Young's modulus using the Sneddon model for forces > 200 pN. We neglect the information at low indentations, since according to [64], the Sneddon model is accurate at higher indentation δ, allowing us to extract the cortical Young's modulus E c as: where F is the measured force. Assuming a Poisson's ratio ν of 0.5 [65], we can fit this formula to the typical force-indentation curves obtained by AFM for every pixel on the cell's surface (we use the Levenberg-Marquard algorithm in MINPACK through its python-interface provided by SciPy for curve-fitting). To extract the stiffness of the cortical layer, we select regions on the cell away from the nucleus where the average cell height is very low so that we can assume that the measured stiffness is indeed the compressive stiffness of the cell's cortex and not dominated by effects from bending of the cortical layer or the intra-cellular fluid-see S1 Text and S2 Fig. To limit the influence of the underlying substrate, the maximal force was chosen to keep the indentation depth to less than 20-30% of the height of cortex cell thickness.
The full procedure to select usable patches within the AFM stiffness maps is detailed in the supplementary information. The global average over all measured cells and all patches yields an estimated cortical stiffness E c = 3.5 ± 2 kPa (Fig 2). Calibration of cell mechanical model For a relatively thin cortical 'sheet' , and assuming that the cortex consists out of some homogenous elastic material, the stretching stiffness k s and bending energy k b can be related to the cortical Young's modulus E c and the cortical thickness t c : where we usually assume the Poisson's ratio of the actin cortex v c to be close to 0.5 [66]. Having determined the effective stiffness of the cortical shell and its thickness from AFM, these equations allow us to calculate the parameters of the mechanical cell model. To evaluate our procedure, these estimated mechanical parameters can be compared to simulated Micropipette Aspiration (MA). Hereto, a simulation was set up where a spherical cell is aspirated into a thin cylindrical structure with a rounded tip and radius R p - Fig 3. The relationship between the applied under-pressure in the pipette and the aspirated length L p of the cell expresses an effective equilibrium Young's Modulus E 1 which can be compared to experimental values obtained using the same technique [67]: where Φ % 2.1 is a scaling factor. The cell's Young's modulus obtained by applying this procedure-see Fig 3-

Design of simulations
The flow profile around a single cell attached to the scaffold is computed by solving the immersed boundary problem at the scale of the investigated cell. For this purpose, the Eulerian computational domain O corresponds to a box of a few hundred microns wide/long containing the cell and not the whole scaffold pore-see Fig 4F.
In order to obtain the magnitude of flow velocity which is to be used as a Dirichlet boundary condition on the microscopic domain O, see Eqs (1) and (2), Stokes' equation was solved on an Poisson's ratio nucleus v n 0.5 - [69] doi:10.1371/journal.pcbi.1005108.t001

Implementation
The Immersed Boundary implementation was realized using the Finite Element software Free-FEM++ [70], which solves the Stokes flow problem, with the Lagrangian forces computed in a coupled module implemented in the particle-based simulation platform Mpacts [71].

Ethics statement
This procedure was approved by the ethics committee for Human Medical Research KU Leuven (ML7861). Patient informed written consent was provided by the legal guardian.

Results and Discussion
Four potentially relevant mechanical measures were computed for cells experiencing flow conditions inside a scaffold pore: nodal displacement, cortical tension, normal pressure and local shear stress. Moreover, we compared between four different geometries as illustrated in Fig 5. All these geometries are regularly encountered in experimental set-ups where cells are attached to titanium scaffold struts-see Fig 5. A parameter study was performed to investigate the effect of changes in volumetric flow rate in the bioreactor on the specified mechanical measures and shows a linear relationship for all four quantities (see supplementary information). In the following section, we will discuss for each of these geometries the effect of flow on the four mechanical measures. In order to permit a direct comparison, the spatial distribution of every quantity will be shown for each geometry using the same scaling. Fig 6 summarizes the effect of flow on the nodal displacements. Very low displacements were obtained for flow parallel to the cylindrical strut (A), while intermediate displacements were found for flow perpendicular to the cylindrical strut (F) and relatively high displacements were observed for the cell forming a bridge at a scaffold corner (B).
For inlet bioreactor flow rates in the order of 1 ml/min, maximal displacements are in the range of 30 nm which has been reported as a critical displacement for the detachment of mesenchymal stem cells of bridged morphology from irregular scaffolds [72,73]. However, these studies were carried out using a one-way fluid-structure interaction (FSI) approach and did not consider the influence of cellular deformation on the surrounding fluid flow, something that in this study was included (two way interaction between cell and fluid). In another study [74], 'microdisplacements by machine vision photogrammetry' (DISMAP) was used to measure flow-induced strain of osteocytes on a flat substrate. Here, a linear strain-shear stress relationship was found. In our simulations, configurations (F) and (A) produce a maximal cell strain of resp. 0.412% and 0.062%, for empty scaffold shear stress magnitudes of resp. 0.0723 Pa and 0.0143 Pa. Using the linear relationship found in [74], strains of resp. 0.307% and 0.0608% would be expected. Even though the cell types are not identical, very similar values would be obtained by our computational model as found by these experiments. Moreover, the higher strain found in the (F) configuration is an expected result of the less shielded cell location.
The observed ranges of cell deformation are much smaller than typical deformations of cells in tissue: e.g. for chondrocytes in mature cartilage, strains higher than 20% (i.e. more than 1 μm) have been measured upon tissue compression [75]. However, it should be stressed that our simulations report the instantaneous elastic response in deformation to a step increase in flow velocity, and neglect the viscous deformations that might occur when cells are exposed to constant flow conditions for a long time (i.e. days).
In Fig 7, the distribution of cortical tension T is shown for the four different geometrical configurations. Positive values of T indicate tensile conditions, whereas negative values of T indicate compressive stresses in the cortical shell. Unlike the cell deformations, which were maximal for the cell bridging between two struts (B), the maximal tension T occurs for cells on cylindrical struts with flow perpendicular to the strut (F) and (T). Moreover, maximal cortical tensions are observed close to the nucleus and close to FAs, with tensile stresses occurring at the side of incoming flow and compressive stresses at the side of out-going flow. Maximal tensions, which are highly localized, are around 5 pN/μm. For comparison, these values are several orders of magnitude below the values for membrane rupture [76]. Other experimental studies have looked at the induction of blebbing, for which it was reported that cortical tensions of at least 200 pN/μm were required [77], while inside blebs, cortical tensions between 10 and 100 pN/μm were measured [78]. The cell's acto-myosin contractility alone creates an average resting cortical tension in the order of 0.5pN/μm [79]. In other words, the predicted additional cortical tension due to shear flow is relatively low, but it cannot be excluded that these tensions could nonetheless result in some conformational changes in the cytoskeleton. In vivo, cells in the bone tissue have been found to experience shear stresses of 0.8-3.0 Pa during routine physical activity [41,80]. A maximal shear stress value of 0.77 Pa was observed on the surface of osteocytes due to flow in the pericellular domain [44].
Compared to the expected wall shear stresses in an empty scaffold, we see that the maximum stress predicted by the IBM is generally about twice as high, while the average stress over the complete cell surface is significantly lower. This clearly indicates the importance of taking the shape and mechanical response of the cells into account for estimating their relevant wall shear stress. Regarding the (F) and (T) configurations, as expected, the distribution of shear stress is more homogenous with half of the cell surface exposed to a higher value than 0.015 Pa, while the two other configurations present most of their surface exposed to low shear stress value (below 0.005 Pa). For the flow rate level used in this study, we have previously reported the effect of increasing flow rates resulting in osteogenic priming of hPDCs in the absence of supplementary growth factors [8] with genes such as osterix and bone sialoprotein being slightly upregulated. Additionally, again for similar flow rates and in the presence of osteoinductive medium, hPDCs have been shown to secrete higher levels of ECM and to enhance mineralization for increasing flow rates [81]. It has been observed [82] that for a shear stress value exceeding a threshold value of 0.088 Pa, human MSCs can detach from the scaffold surface and by doing so negatively affect the final properties of tissue engineered constructs, resulting in an inhomogeneous distribution of neotissue across the scaffold [83]. This illustrates the importance of having a numerical model such as the one presented in this study, to quantify and characterize the mechanical environment that cells experience in novel scaffold designs in order to avoid the development of suboptimal or detrimental mechanical regimes.

Conclusions and outlook
In the presented work, a novel application of the immersed boundary method was developed, representing a deformable cell exposed to microscopic flow and attached to a 3D scaffold inside a perfusion bioreactor. Cells were represented by a deformable Lagrangian surface mesh, which was immersed in an Eulerian fluid domain, with flow in the Stokes regime. We demonstrated the effect of shear flow for multiple realistic geometrical cell configurations and strut locations inside a regular pore scaffold. This tool can be used to estimate shear flow conditions directly on the surface of individual cells, and assess the micro-scale variability of mechanical conditions inside single scaffold pores. The instantaneous cell stiffness was measured using AFM experiments on hPDCs, and the mechanical model was calibrated using micro-pipette aspiration simulations in the range of short term deformations. Simulations confirmed that mechanical cues originating from the flow are highly dependent on the exact geometry of the cell and its environment. For example, a cell on a cylindrical strut with flow perpendicular to the strut will experience a much larger shear stress than a cell on a similar strut with parallel flow. This should be a major consideration when designing novel scaffold designs. Moreover, it was found that wall shear stress calculated in the empty scaffold would underestimate the actual maximal wall shear stress experienced by the cells by a factor of two in the investigated cases. Next, the model was used to estimate the additional instantaneous flow-induced cell deformation, tension and pressure. Compared to the cell-generated deformation and tension due to acto-myosin activity these values are very small for the applied realistic flow conditions. For example, the predicted magnitude of additional cortical tension due to flow is much lower than the pre-stress of a spread out cell. Hence, it is not expected to affect the tensional homeostasis of a cell in steady-state conditions. Nonetheless, in cyclic conditions, small deviations from the resting state can still have pronounced biological effects, as is explained by the tensegrity model [32,33]. In the same vein, small cytoskeletal deformations away from the cell's resting configuration might alter the configuration of several intracellular mechanosensing molecules, which are linked to downstream targets in pathways such as the mitogen-activated protein kinase (MAPK) or phosphoinositide-3-kinase (PI3K) pathways [84]. Since the predicted strain of the cell is highest near its attachment points on the substrate, which also constitute focal sites of mechanotransduction machinery, a discernible biological response could result even from small deformations, as long as the dynamics of the perturbations in flow conditions are faster than the relaxation time of tensional homeostasis.
In the current work, a strongly simplified model for the mechanical behavior of single cells was used, which limits its predictive value to small deviations from the resting state. A more indepth computational analysis of the mechanisms at play involved in large cell deformations would require a more detailed model of the mechanoadaptive behavior of the cytoskeleton and will be reserved for future research. Similarly, the effect of shear flow on the long timescale viscous-like deformation of living cells remains to be investigated. This would also require a more elaborate description of the cell's mechanical behavior, which for this study was greatly simplified and limited to linearly elastic deformation. The model's limitation of small deformations and short timescales will mainly restrict in scope the predictions of cell deformation and cortical tension, whereas predictions of shear stress and local pressure, being surface properties, are expected to remain largely unaffected. Finally, to simulate adhesion and detachment behavior (e.g. in very high shear flows), the presented methodology has to be extended since adhesion is only implicitly captured by placing FAs out of the fluid domain thereby fixing them in space independently of applied forces. For this, an adhesion force formulation as proposed in [62] could be included. A parameter study (see supplementary information) showed linear behavior in the relevant cell-mechanical and flow parameters, showing that the model can be used to inter-/extrapolate to different cell types and flow conditions. For a cell of thickness 5 μm facing flow on a strut of diameter 200 μm, the wall shear stress can be estimated as: τ % 0.08ÁQ in . Evidently, the maximal wall shear stress experienced by a cell does not depend strongly on the cell's mechanics. This implies that even though a cell may undergo structural changes (e.g. migration, re-alignment), it can still reliably 'sense' the shear flow (e.g. with its primary cilium). This study constitutes an important step towards model-based control of a cell's biophysical microenvironment (stem cell niche engineering) in a perfusion bioreactor. From Q in the Dirichlet boundary conditions in the micro-scale model were determined using a CFD simulation of the complete scaffold pore- Fig 4A. The resulting maximal deformation, pressure, shear stress and cortical tension were quantified. One might notice that the dependence on Q in is linear, which is due to the Stokes' flow regime, which is valid for the investigated range of flow rates. Except for the maximal deformations, the effect of the cells' stiffness is very small.