The calculation of fluid – induced WSS within empty scaffolds is based on the assumption that the WSS at the scaffold surface is a good representation of the WSS sensed by the cells that are attached to the scaffold surfaces. It also assumes that the cells attach mostly flat to the scaffold surface in the initial state, with a minimal cell volume with respect to the pore volume. This assumption has been shown to be met for some experiments/scaffold materials (Figure 1A), but not for all (Figure 1B). For calculating the mechanical strain in empty scaffolds, it is assumed that the cells are subjected to the strain magnitude at the location of the scaffold that they are attached to (Olivares et al., 2009; Laurent et al., 2014). Such assumptions can be reasonable if the scaffold material is much stiffer than the cell and no substantial ECM has been formed yet.

Porosity is the main determinant for scaffold permeability and thus the amount of flow through the scaffold in perfusion/stirring bioreactors and was found to be an important parameter determining the results of BTE (e.g. seeded cell density, cell proliferation, ECM production, etc.) (Grayson et al., 2008; Panseri et al., 2021). The effect of porosity on the permeability of the scaffold, which describes the amount of flow through the scaffold, can be calculated by the Kozeny-Carman Equation (Eq. 1) (Van Bael et al., 2012; Egan, 2019): κ = ( 1 C k ) ⋅ φ S s 2 (1) where, κ is the permeability, φ is the porosity, c _k is the Kozeny constant and S _s is the specific surface area calculated as the surface area divided by the total volume of the struts.

This equation demonstrates that scaffold permeability linearly increases with porosity. This has also been demonstrated by experimental measurements (Zhang et al., 2019). The relationship between the permeability, the fluid velocity and the WSS, however, is complex and also depends on the pore geometry. Ali and Sen (2018) employed a CFD approach to investigate the influence of porosity on the permeability and WSS, and they found that under a fluid velocity of 0.7 mm/s, the average WSS decreased from 131 to 27 mPa with an increase in porosity from 65 to 90% for the gyroid pore shape (Figure 2H). This trend also happened for a diamond pore shape (Figure 2F) (Ali and Sen, 2018). Melchels et al. (2011) designed and manufactured a scaffold (gyroid pores in Figure 2H) with different porosities (40–85%) in different regions. This resulted in different shear rates (SRs) in the regions with different porosities under perfusion flow, e.g. 10–40 s⁻¹, i.e. higher SR in the region with higher porosity (Melchels et al., 2011).

Porosity is also the main determinant for the amount of strain in the scaffold when an external force is applied according to both micro-FE analyses and experimental characterisation (Hannink and Arts, 2011; Castro et al., 2020). If under pressure/compressive force loading through bioreactors, a scaffold with a lower porosity (i.e. higher overall structural stiffness) will show less deformation, thus cells that attach on the struts will receive lower strain. Not only for overall structural stiffness, the porosity also can influence the local stress/strain concentration under compression loading. For instance, in the computational study by Hendrikson et al. (2017), scaffolds that had a cubic pore shape, a pore dimension of 151 μm, but different porosities (74 vs 64%) were compared in terms of octahedral shear strain under a compressive strain of 10%. It was found that the scaffold with higher porosity (74%) had a larger strut area in the low strain range (e.g. <10%) than the one with lower porosity (64%) (Hendrikson et al., 2017). Similarly to section 2.2 and 2.3, under dynamic compression, the porosity will also have an influence on the resultant WSS that is induced by fluid flow due to compression. Zhao et al. (2016) compared scaffolds with porosities of 60–90% and found that a higher porosity resulted in lower WSS under dynamic compression. By increasing the pore dimension, the influence of porosity on WSS became smaller, e.g. for cubic pore, dimension = 100 µm: average WSS = 3.5 mPa when porosity = 60%, average WSS = 2.5 mPa when porosity = 90%; for cubic pore, dimension = 300 µm: average WSS = 1.5 mPa when porosity = 60%, average WSS = 1.1 mPa when porosity = 90% under dynamic compression (strain = 1%, frequency = 1 Hz) (Zhao et al., 2016). So, under dynamic compression, the lower porosity and pore dimension can result in higher fluid flow-induced WSS.

Pore dimension is the main factor that determines fluid-induced WSS under perfusion flow (Fu et al., 2021). Also, pore dimension is one of the factors that can influence cell attachment (e.g. flatly attached on struts/bridging over struts in Figure 1) (Guo et al., 2015; Yamashita et al., 2016). Previous mechanobiological studies have investigated cell responses via tuning the scaffold pore dimensions while keeping the pore shape constant (Bartnikowski et al., 2014; Ouyang et al., 2019). For example, Bartnikowski et al. (2014) quantified the WSS within the scaffolds that had cubic pores (with rounded profile in Figure 2G) and a porosity of 60% but with different pore dimensions (625 vs 1250 µm). It was found that the scaffold with smaller pore dimension provided a higher WSS: maximum WSS = 1979 mPa/average WSS = 500 mPa (pore dimension = 625 µm) vs maximum WSS = 837 mPa/average WSS = 120 mPa (pore dimension = 1250 µm) under a flow rate of 1 ml/h (0.61 µm/s). These scaffolds were then applied in an in vitro cell experiment where it was found that the DNA amount was significantly higher in the cell seeded scaffold with larger pore dimension (1250 µm) (Bartnikowski et al., 2014). Whereas the accurate calculation of the actual WSS requires performing a CFD analysis for the (often complex) scaffold pore geometries, simple mathematical equation can be used to estimate the WSS (Zhao et al., 2016): τ a v = a 1 ⋅ ( d L c ) b 1 (2) Where, τ _a is the average WSS within the scaffold, v the applied fluid velocity, d the pore dimension, L _c is the characteristic length (L _c = 1 µm), and a ₁ and b ₁ constants that depend on the pore shape and porosity. A limitation of this approach is that the equations are only a good approximation for a limited set of pore shapes (i.e. cube with flat profile and sphere in Figures 2C,E) and porosities (i.e. 60%–90%).

Dynamic compression not only generates mechanical strains in the struts but can also result in WSS on the strut surfaces, which was mostly ignored in previous mechanobiological studies regardless of the compressive strain magnitude or frequency. It was found that the resultant average WSS was proportional to the applied compressive strain (Milan et al., 2009). Moreover, it was highly dependent on the scaffold pore dimensions (Zhao et al., 2016). Here also, simple equations were introduced to estimate the WSS due to compression (Zhao et al. (2016)): τ a ε a p p = a 2 ⋅ ( d L c ) b 2 (3) where, τ _a is the average WSS within the scaffold, ε _app is the applied compressive strain by bioreactor, d the pore dimension, L _c is the characteristic length (L _c = 1 µm), and a ₂ and b ₂ constants that depend on the pore shape and porosity. Similar as in Eq. 2, one of the limitations of this approach is that the equations are a good approximation only for a limited set of pore shapes (i.e. cube with flat profile and sphere in Figures 2C,E) and porosities (i.e. 60%–90%). Also, this Eq. 3 is only applicable for a dynamic compression frequency of 1.0 Hz and needs to be adapted to other frequencies. For all other cases beyond the aforementioned ones, a FSI analyses will be needed to accurately calculate the WSS.

Under mechanical compression, the stress/strain distribution can be influenced by the pore dimension. Ribeiro et al. (2017) investigated this based on the scaffolds with pore dimensions of 740 and 370 µm using an FE model. There, an unconfined compression loading with a strain of 15% was applied on both scaffolds. Their results showed that the maximum value of compressive stress was similar between the two scaffolds with different pore dimensions, i.e. maximum compressive stress = 27.7 MPa in pore dimension of 740 µm vs maximum compressive stress = 25.9 MPa in pore dimension of 370 µm. However, the scaffold with larger pores (pore dimension = 740 µm) had more regions (area) with higher stress than the one with smaller pores (pore dimension = 370 µm) (Panadero et al., 2015).

The effect of pore shape on the fluid-induced WSS in flow perfusion and spinner flask bioreactors is difficult to predict. Some scaffolds with different regular pore shapes (but same pore size and porosity) have similar WSS, but some do not. The commonly (designed) regular pore shapes include sphere, cube, gyroid, prism, etc. (Figures 2A–H), which can be manufactured by 3D printing/additive manufacturing technology. Prochor and Gryko (2021) quantified the WSS within scaffolds that have different regular pore shapes (e.g. triangular prism with rounded and flat profiles, cube, octagonal prism and sphere in Figures 2A–E) under perfusion flow. It was found that the scaffold with triangular prism (with rounded profile) experienced the highest WSS, whereas the scaffold with spherical pores experienced the lowest WSS under the same flow rate. The maximum WSS within spherical pores and cubical pores were identical. However, this can be different within different bioreactors that generate fluid – induced WSS. In a combined experimental and computational BTE study by Rubert et al. (2020), the average WSS within the scaffolds with spherical pores (average pore diameter = 330 μm, porosity = 84.7%) and cubical pores (average pore diameter = 330 μm, porosity = 92.8%) were 0.42 and 0.81 mPa respectively in a spinner flask (70 RPM). This was associated to upregulated osteoblast cell differentiation and ECM formation within cubic pores, while ECM mineralisation was enhanced within the spherical pores (Rubert et al., 2020).

Porous scaffolds also can have irregular pore shapes, which are for example obtained from more traditional fabrication methods such as porogen leaching (Figure 2I). Studies have found that the irregularity of the pore shape does not have a distinct influence on the fluid – induced WSS, once the pore dimensions and porosity are similar. For example, Zermatten et al. (2014) investigated the influence of the pore irregularity on the internal WSS using scaffolds with regular cubical pores (with rounded profile) and highly irregular pores (Figures 2G,I). Although the other two parameters, pore dimensions (regular: 0.22 mm vs irregular: 0.16 mm) and porosity (regular: 38% vs irregular: 55%) were not exactly the same, the average WSS within irregular and irregular pores have high similarity (regular: 3.08 mPa vs irregular: 3.68 mPa) under a perfusion fluid velocity of 0.066 mm/s (Zermatten et al., 2014). One limitation of simulating the micro-fluidic environment within these highly irregular pores (at whole scaffold level) was the high computational cost (Santamaría et al., 2013; Zermatten et al., 2014). To address this challenge, Zhao et al. (2019) developed a more versatile technique, creating a multiscale and multiphasic CFD model. In this approach, small but representative parts of the scaffold are being used for the generation of a microstructural model of the pore environment, which are then coupled with a macro-model representing the whole scaffold in which the microstructure is homogenised. The macro-model can be used to calculate the fluid flow at larger length scales, that then can be applied to the micro-model to calculate the local WSS at the cell level. As only small parts of the scaffold need to be modelled in detail, this approach can reduce the computational costs while still providing results at the cell-level. It has been shown that with this multiscale and multiphasic CFD model, calculations of resulting WSS within any scaffold with highly irregular pore shape is possible even using a normal computer (e.g. 16 GB RAM, Intel i7 CPU). However, this technique requires that the Reynolds (Re) number should less than 1 when using Darcy’s law for homogenisation (Chaudhary et al., 2011).

Under mechanical compression, the pore shape can have a distinct influence on the overall structural stiffness of the scaffold (Castro et al., 2020; Jahir-Hussain et al., 2021). According to their calculations, the triangular pore shape resulted in the highest structural stiffness of the scaffold and the spherical pore shape resulted in the lowest stiffness among the various pore shapes (spherical, cubic, hexagonal and triangular). Under pressure/compressive force, scaffolds with a lower structural stiffness (e.g. with spherical pores) will have a higher strain in the struts than the ones with a higher structural stiffness (e.g. with cubic/hexagonal/triangular pores). This difference will translate to differences in strain sensed by cells attached to the struts. In cell culture experiments applying compression to stimulate cells, usually dynamic compression is applied. As mentioned in sections 2.2 and 2.3, this dynamic compression also generates WSS within the pores of the scaffold. A FSI approach for quantifying the WSS during dynamic compression has found that the WSS was higher within spherical pores than that within cubical pores, e.g. 5.5 mPa within spherical pores (pore diameter = 100 μm, porosity = 60%) and 3.5 mPa within cubical pores (pore diameter = 100 μm, porosity = 60%) under an applied compressive strain of 1.0% and at a frequency of 1.0 Hz (Zhao et al., 2016). Therefore, to precisely quantify the WSS due to dynamic compression, the pore shape needs to be explicitly reflected in the computational model.

Some scaffolds have extremely anisotropic pores, such as those with unidirectional channels or holes, as shown in (Deville et al., 2006, 2007; Munch et al., 2009; Pourhaghgouy et al., 2016). For this type of scaffolds, the above discussed influence of porosity, pore dimension and pore shape on the internal micro-mechanical environment is still applicable for external loading in the unidirectional orientation (e.g. fluid perfusion/unidirectional mechanical compression/stretching), but not for external loading in non-unidirectional directions (e.g. spinner flask/non-unidirectional compression/stretching).

In some BTE experiments, dynamic cell seeding is used for improving the seeding efficiency and/or distribution of seeded cell in the porous scaffolds. Perfusion flow is usually used for dynamic seeding. During this process, the fluid force can also mechanically stimulate the cells through cell deformation which consequently can promote cellular processes (Rüberg and Aznar, 2016; Serrano-Alcalde et al., 2017).

For seeded cells, Jungreuthmayer et al. (2009) and Mccoy et al. (2012) modelled cells as flatly attached and as bridged morphologies within collagen glycosaminoglycan (CG) scaffolds, which had irregular pore shapes (Figure 3A). It was found that the influence of cell morphology (attached/bridged) on the cellular WSS depends on the locations within scaffolds (Guyot et al., 2016b). Furthermore, in the study by Mccoy et al. (2012), three CG scaffolds with different pore dimensions (85, 120 and 325 µm) but equal porosity (90%) were compared in terms of resultant WSS on cells. It was found that the average WSS on all cells (both bridged and attached morphologies) was 165, 176 and 155 mPa, respectively for the pore dimensions of 85, 120 and 325 µm under a perfusion fluid velocity of 235 µm/s, and the WSS was proportional to the fluid velocity (Mccoy et al., 2012).

To investigate the influence of scaffold pore geometry on the changing micro-mechanical environment over time, tissue growth models have been introduced. These are coupled with FE/CFD models. To account for tissue growth, various studies have developed mathematical models to describe neo-tissue formation assuming appositional growth in static conditions (Nava et al., 2013; Guyot et al., 2014; Herklotz et al., 2015; Egan et al., 2018). In the mathematical model, the cells and ECM are usually homogenised (Figure 3B). One of the most commonly used models for describing the tissue growth kinetics is based on the level set (LS) method, which is available in both commercial and open-source software packages (e.g. COMSOL, ANSYS, OPENFOAM, etc.). It allows to model appositional tissue growth as illustrated in Figure 3B. The governing equation is (Guyot et al., 2014): ∂ ψ ∂ t + ( V G ⋅ n Γ ) ⋅ ∇ ψ = 0 (4) where, n _Γ is the normal unit vector to the interface between neo-tissue and medium domains, ψ is the LS function and V _G is neo-tissue growth velocity governed by local the struts curvature κ _c in Eq. 5: V G = { − λ ⋅ κ c           ( κ c > 0 ) 0                               ( κ c ≤ 0 ) (5) Where, λ is the tissue growth rate.

To investigate the influence of the micro-mechanical environment changes during neo-tissue growth within scaffolds that have different pore geometries, the tissue growth model needs to be coupled with the FE/CFD model by introducing a WSS-dependent term into Eq. 5. The WSS (τ) in Eqs 6, 7 by Guyot et al. (2015, 2016a) then is computed by a CFD model:   V G = {   − A ⋅ κ ⋅ f ( τ )           κ   >   0 0                                                 κ   ≤   0 (6) f ( τ ) = { 0.5 + 0.5 τ a 1     0 ≤ τ ≤ a 1 1                                       a 1 ≤ τ ≤ a 2 τ − a 3     a 2 −     a 3       a 2 ≤ τ ≤ a 3         0                                   τ ≥ a 3 (7) Where, a ₁ and a ₂ are the minimal and maximal shear stresses enhancing neo-tissue formation and a ₃ the critical shear stress.

Then the computational model is applied to scaffolds that have different pore geometries. For instance, Guyot et al. (2015) applied the model on scaffolds with two different pore shapes (i.e. cubic shape with rounded profile and pore dimensions of 650 µm vs diamond shape and pore dimension of 750 µm in Figures 2F,G). It was found that under the same amount of tissue produced within the scaffolds higher shear stress occurred in the neo-tissue within the scaffold with diamond shape than that in the scaffold with cubic shape, e.g. when 30% of the porous volume was filled with neo-tissue, the average shear stress in the neo-tissue was 175 mPa within the diamond pores, while it was 80 mPa within the cubic pores (Guyot et al., 2015).

One of the limitations of these computational models is the uncertainty of the parameter values (such as λ, a ₁, a ₂ and a ₃ in Eqs 5–7). As these are empirically determined constants, they may need to be changed depending on parameters that influence tissue formation (e.g. the number of cells in the culture, the type of cells, scaffold-related attachment of cells, the density of the deposited ECM and whether or not it is mineralised). Whereas after fitting these constants to experimental results these equations thus may well describe the effect of changes made within that specific experiment, they may not well describe the outcome of other experiments. To reduce the number of parameters in the tissue growth model, recent computational studies have employed second order diffusion equations to model tissue growth kinetics (Buenzli et al., 2020; Zhao et al., 2020a). The main advantage of using this diffusion equations over the LS method is that fewer parameters need to be determined. For example, diffusion equations can already model the curvature – dependent tissue growth without adding the curvature parameter κ in the equation as that in LS method (Buenzli et al., 2020). Therefore, in modelling the scaffold pore geometry for tissue growth kinetics, if the curvature is not a parameter that needs to be explicitly assessed, a computational model based on a diffusion equation will be a good choice. Otherwise, a computational model based on LS method is suggested. Another limitation is that these computational models assume appositional tissue formation towards the centre of the pores. In reality, however, interstitial formation, in which the tissue is infiltrated within the pores rather than being attached on the struts surfaces is also observed in many cases (Li et al., 2009) (as illustrated in Figure 3C). The resultant WSS on cells under interstitial tissue formation was quantified and compared to appositional tissue formation (Zhao et al., 2020b). Distinct difference in WSS between two cases were found, even if the same amount of newly formed tissue was present. This implies that computational models that assume appositional tissue growth cannot well predict the micro-mechanical environment in case of substantial interstitial tissue formation. Estimating the influence of scaffold pore geometry on the micro-mechanical environment by taking the tissue into account also needs to consider whether the cell/tissue growth is appositional or interstitial. Even then, this remains challenging due to the high variability in tissue formation.

Different from LS method and diffusion equation, some other computational studies employed a simple voxel – FE based method to simulate the tissue growth within scaffolds (Adachi et al., 2006; Nasello et al., 2021). In this method, modelling the neo-tissue generation within scaffolds was achieved by adding elements on the scaffold surfaces according to the applied stress in elements where the cells are located. Therefore, this voxel – FE based method does not require mathematical functions for tissue growth kinetics. However, this method is limited to simulate the neo-tissue growth under mechanical stimulation only, and not under static conditions.