Page 290 - IJB-10-3

P. 290

International Journal of Bioprinting Design and optimization of 3DP bioscaffolds

time. This model provides an intuitive description of the 2.5. Numerical method
relationship between cell density changes over time and A single timestep solution procedure for the multi-physics
the cell growth rate. Such a step-wise profile is not only model is illustrated in Figure 2. Specifically, the parameters,
suitable for describing the logarithmic phase of cell growth such as flow rate (velocity), oxygen diffusion coefficient,
but also can effectively incorporate the phase lag. initial cell density, initial oxygen concentration and initial
porosity, were fed into the model. Computational fluid
Furthermore, to take into account the variation in
the scaffold’s porosity resulting from the increase in cell dynamics (CFD) simulations for the laminar and porous
media flow were then performed to predict the oxygen
number, the following relation is coupled to the model:
mass transfer in the system where oxygen convection and
diffusion occurred. Subsequently, the obtained oxygen
ε()t = ε − V ⋅ ρ ()t (22) concentration results, such as the oxygen concentration,
0
c
c
cell density, and other parameters of each node in the
where ε represents the initial porosity of the medium;
0
V denotes the volume of a single cell; and ρ is the number structure, as the inputs, were fed to the cell oxygen
consumption model to calculate the oxygen consumption,
c
c
of cells per unit volume. yielding new oxygen concentration values for simulating
2.4. Initial values and boundary conditions cell proliferation. These steps were repeated to obtain
The initial values of oxygen concentration (C ), cell continuous results for the parameters within the target
0
density (ρ ), scaffold porosity (ε ), and oxygen diffusivity time interval.
0
0
(D ) in the system and model boundary conditions are To solve the multi-physics model more efficiently and
0
summarized in Table 1. Boundary 1 is the inlet of the conveniently, COMSOL Multi-physics software was used
nutrient solution flow, as depicted in Figure 1a. Both the in this study. Given that various physical fields, such as fluid
initial flow rate of the fluid and oxygen concentration dynamics, mass transport, and cell kinetics, were coupled,
at boundary 1 were chosen to be constant, being appropriate finite element mesh generation should be
Q = 600 µL/min and C = 222.5 µM, respectively. obtained to balance solution accuracy and computational
0
in
The velocities and oxygen concentrations should be efficiency. Due to the different physical processes involved
continuous at the interfaces (boundaries 3, 4, and 5) in each domain, a partitioned mesh approach was adopted
between the scaffold and nutrient solution flow, as for mesh generation, as illustrated in Figure 3a. Region 3
indicated by the relation in Table 1. At boundary 2 indicates the porous media flow domain in the channels
where the flow is in contact with the chamber wall, the of the biological scaffold, which is coupled additionally
boundary was set to be a no-slip condition, while the with oxygen consumption and cell growth processes.
oxygen concentration was set to be a no-flux condition. Hence, the mesh density and refinement in Region 3 must
At outlet boundary 6, the fluid satisfies a zero gauge be higher to resolve the underlying physics adequately,
pressure condition. while the mesh in regions 1, 2, and 4, where laminar flow

Table 1. Initial values and boundary conditions for model set-up

Initial values Boundary conditions

Volume 10 Issue 3 (2024) 282 doi: 10.36922/ijb.1838

285 286 287 288 289 290 291 292 293 294 295