Mathematical Modeling Reveals the Role of Hypoxia in the Promotion of Human Mesenchymal Stem Cell Long-Term Expansion

Many experimental studies have found that human mesenchymal stem cells (MSCs) in long-term culture exhibited enhanced cell proliferation and prolonged lifespan under hypoxia (around 1%–7% oxygen) against the normoxic condition (about 21% oxygen). Inspired by the experimental findings, we aimed to investigate the hypoxic effects on MSC expansion quantitatively through mathematical modeling to elucidate the corresponding biological mechanism. A two-compartment model based on ordinary differential equations (ODEs), which incorporate cellular division and senescence via state transition, was developed to describe the MSC expansion process. Parameters of this model were fitted to experimental data and used to interpret the different proliferative capacities of MSCs under hypoxia and normoxia along with model sensitivity analysis. The proposed model was tested on data from two separate experimental studies, and it could reproduce the observed growth characteristics in both conditions. Overall, this compartmental model with a logistic state transition rate was sufficient to explain the experimental findings and highlighted the promotive role of hypoxia in MSC proliferation. This in silico study suggests that hypoxia can enhance MSC long-term expansion mainly by delaying replicative senescence, which is indicated by the slowdown of the state transition rate in our model. Therefore, this explanatory model may provide theoretical proof for the experimentally observed MSC growth superiority under hypoxia and has the potential to further optimize MSC culture protocols for regenerative medicine applications.

can be found in Table S1, of which the midpoint time is annotated in the figure as and ℎ for normoxia and hypoxia respectively.  Table 1 as a contrast of Figure 4 in the main text. As shown here, even if we assume the net expansion rate 1 is slightly larger in the normoxic condition, there is a still a large gap between the final population size under the two conditions.
Besides, the fitting results change quite little in comparison with Figure 4. Thus, we can conclude that the minute difference of the fitted 1 value under two oxygen conditions plays just a negligible role in causing the significant difference of MSC expansion efficiency.

Implementation details of nonlinear regression in MATLAB
In this section, we present more details about how to solve the nonlinear regression problem (5) corresponding to Figure 3 in the main text, i.e., how to minimize ( ). Once we choose an initial value 0 for the parameter vector, we can simply feed the data into mature optimization algorithms in MATLAB, such as fminsearch, fmincon or lsqcurfit. However, it is well known that the optimization of a nonlinear, nonconvex objective function may be sensitive to initial values due to the existence of multiple local optima. Therefore, in the following we mainly explain how to choose a good initial value for the parameter vector . In our study, since for each experiment there are only few data points, we simply use an exhaustive grid search method to determine the initial parameter values [1], detailed as follows.
Step 1: Determine a minimal parameter space, according to our prior knowledge or by observing the experimental data.
In equation (5) of the main text, the constraints for the parameters are derived from merely mathematical requirements. However, with some prior biological knowledge or simply by observing the experimental data, we can further narrow the parameter space.
For example, with the data points of Experiment A presented in Figure 4, it is easy to estimate the slope of the first stage, i.e., the exponential growth stage, to be around 0.5. According to equation (8), we know that the theoretical slope is 1 log 2 and thereby the value of 1 should be less than 0.5, since log 2 ≅ 1.44 > 1. As for the parameter 20 , which represents the death rate, it should be a small value according to our biological knowledge of MSC proliferation: we simply guess that 0 < 20 < 0.2. It is also easy to determine a rough range for the parameter , the midpoint time of the logistic function, by noticing that in Figure 4 the exponential growth doesn't stop until day 50 and the plateau stage is reached around day 80 for both conditions. Thus, we know approximately 50 < < 80. For the other two parameters, and , it is a little difficult to further narrow their range, and we preserve their upper bound listed in (5), i.e., 1 < < 1 and 0 < < 1.
With such a narrowed parameter space, it is easier to choose sensible starting-values for the optimization.
Step 2: With the parameter space determined in step 1, perform a coarse-grained grid search to try different initial values for the optimization routine.
The classical grid search method constructs a discrete grid based on the parameter space composed of all parameters. Usually, the points are sampled equidistantly as the initial values to be tested for optimization. In our implementation, for this 1 st -round coarse grid search, we pick 5 equidistant values in 8 the range of each parameter. Therefore, for the total 5 parameters, there are 5 5 = 3125 trials made in total. That is, we have tested 3125 initial parameter values uniformly sampled from the parameter space, and for each initial value the optimization routine is executed and yields a minimizer ̂= min ( ), which depends on the initial values 0 . Define the best minimizer in this round as the one leading to the minimum objective value: ̂ * = min (̂).
Step 3: Perform a fine-grained grid search around the best minimizer we obtained in step 2.
Since the number of grid points increases quickly with model dimension, in the above step 2 we only perform a coarse grid search with 5 samples for each parameter, and a further fine-grained grid search is made in this step. After we get the best minimizer ̂ * in step 2, we build a new grid with the neighboring points of the previously optimal solution ̂ * . Suppose the optimal solution of a parameter obtained in step 2 is * and the original range of this parameter is , then a reduced range is specified to be = [ * − 10 , * + 10 ], which should of course respect the fixed constraint in equation (5). This range reduction principle applies to all the parameters. Afterwards, another grid search routine like the one in step 2 is performed again in this reduced parameter space to further refine the optimal solutions.
The purpose of this fine-grained grid search step is to make it more likely that the global optimum is found.
Even not, after the grid searches in step 2 and step 3, we have much confidence that a local minimum close enough to the global one can be obtained. Finally, the result we report is the best local minimizer found so far, * , which leads to the minimum objective value among all trials in both steps.
In step 2 & 3, given an initial parameter vector, we used the nonlinear programming solver fminsearch in MATLAB to solve the nonlinear optimization problem. Particularly, the most important options are MaxIter=1000 (maximum number of iterations), TolFun=1e-4 (termination tolerance on the function value change) and TolX=1e-4 (termination tolerance on independent variable step). The optimization routine was actually finished in less than 700 iterations for most initial values. We have also tested other optimization options and other solvers such as fmincon. The obtained optimal solutions remains almost the same.

Remarks
1) It is known that the number of grid points increases exponentially with the dimension (number of parameters) of the optimization problem [1]. However, in our study, there are only five 9 parameters and the size of the dataset in each experiment is also very small. Thus, it is still practical to adopt the grid search method to attempt a good starting point for optimization.
2) To repeatedly run the optimization solver with different initial values (starting points), we can use a for-loop manually or resort to the MultiStart algorithm in MATLAB to execute these iterations automatically. Please refer to the online documentation for more details. https://www.mathworks.com/help/gads/multistart.html.
3) There is no technical reason to limit the solvers to be fminsearch or fmincon. Other solvers capable of nonlinear programming can also be used. Furthermore, global optimization solvers like genetic algorithm, differential evolution and particle swarm solvers are also good alternatives, though they may need a higher cost of computation. Interested readers may refer to [2] and [3].