Capturing Mesenchymal Stem Cell Heterogeneity during Osteogenic

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Capturing Mesenchymal Stem Cell Heterogeneity during Osteogenic Differentiation: An Experimental−Modeling Approach Romuald Győ rgy,†,∇ Michail E. Klontzas,‡,§,∥ Margaritis Kostoglou,⊥ Nicki Panoskaltsis,§,# Athanasios Mantalaris,*,‡,∥ and Michael C. Georgiadis*,† †

Department of Chemical Engineering, Aristotle University of Thessaloniki, Thessaloniki, 54124, Greece Biological Systems Engineering Laboratory, Department of Chemical Engineering, Imperial College London, London SW7 2AZ, United Kingdom § School of Medicine, Emory University, Atlanta, Georgia 30332, United States ∥ Wallace H. Coulter Department of Biomedical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, United States ⊥ Department of Chemistry, Aristotle University of Thessaloniki, Thessaloniki 54124, Greece # Department of Haematology, Imperial College London, London SW7 2AZ, United Kingdom

Downloaded via UNIV OF SOUTHERN INDIANA on July 30, 2019 at 01:26:14 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: Currently, there is an increasing clinical demand for bone implants for several orthopedic and maxillofacial surgical procedures. Meeting these demands can be achieved with stem cells differentiated in vitro toward the osteogenic lineage to produce cells suitable for implantation in critical size bone defects. Toward this end, we present a population balance model (PBM) of the osteogenic differentiation of umbilical cord mesenchymal stem cells. Our experimental−modeling approach describes osteogenesis with a dynamic cell PBM incorporating the expression of key osteogenic genes and the activity of metabolic pathways involved in osteogenesis. In addition to accurately capturing the phenotypic heterogeneity of osteogenic differentiation, the model also predicts a delay for the onset of differentiation when the cells grow faster, thus a trade-off between proliferation and differentiation over the course of in vitro osteoinduction. The mathematical model from this work sets a solid foundation for future model-based culture or reactor optimization.

1. BACKGROUND

banks, superior proliferative capacity compared to adult MSCs, and enhanced immunomodulatory properties.5 During in vitro osteogenic differentiation, MSC cultures are heterogeneous cell populations consisting of cells at various stages of the differentiation process.6 Current analyses of differentiating cultures (including gene expression analysis, extracellular matrix production analysis, metabolic analysis, etc.) assume culture homogeneity, disregarding the heterogeneous nature of MSC cultures. Therefore, further understanding of how cell populations evolve over the course of osteogenic differentiation will enable culture optimization and the production of high-quality bone grafts. Osteogenic differentiation is characterized by important phenotypic changes, which are necessary to acquire the osteoblastic physiology. The transition from glycolytic metabolism to energy production through oxidative phosphorylation has been proven to be one of the characteristic changes in MSC

Bone loss during surgical procedures such as major trauma repair, revision joint replacement, tumor resection, congenital skeletal abnormality correction, dental implant fixation, and osteomyelitis treatment all require the use of bone grafts. Autografts are currently used in clinical practice to cover the demand for bone grafts, despite drawbacks such as donor site morbidity, risk of infection, and limited availability. Tissue engineering promises a viable alternative to the use of autografts by combining biomaterials with stem cells and osteoinductive cues for off-the-shelf use of bone grafts in the clinical setting.1,2 The production of high-quality engineered grafts requires efficient stem cell differentiation toward the osteogenic lineage, which guarantees the secretion of extracellular matrix proteins and the successful integration of the implanted graft with the surrounding tissues.3,4 Mesenchymal stem cells (MSCs) have been widely utilized for the production of tissue-engineered grafts, because of their capacity to differentiate into osteoblasts and the fact that they can be isolated from several adult and fetal tissues. Umbilical cord blood (UCB) is an emerging source of MSCs, because of their off-the-shelf availability in cord blood © XXXX American Chemical Society

Received: Revised: Accepted: Published: A

April 11, 2019 July 9, 2019 July 9, 2019 July 9, 2019 DOI: 10.1021/acs.iecr.9b01988 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research physiology over the course of osteogenesis.7−9 Metabolism of UCB MSCs has been found to increase globally during the first 7 days of differentiation and subsequently settle to a more oxidative state that is similar to the metabolism of osteoblasts.9 In order to acquire the osteoblastic phenotype, expression of genes related to differentiation happens in a temporal manner. Such genes include Runt related transcription factor-2 (Runx2), which is expressed early during the differentiation process, and genes related to bone extracellular matrix production, such as SPARC/osteonectin, osteocalcin, and bone sialoprotein, which are expressed at later stages of differentiation.10 Finally, transitioning from a proliferative undifferentiated to a differentiated phenotype is linked to cell cycle changes with a marked gradual increase in average cell cycle duration.11 Stem cell proliferation has been described mathematically in numerous papers using ordinary differential equations (ODEs). Vozzi et al. employed ODEs to investigate the metabolic capacity, senescence hallmarks, and osteoblast-specific gene expression of periosteum-derived progenitor cells, which represent a type of MSC.12 Thalheim et al. utilized a threedimensional (3D) mathematical model of mouse small intestinal crypts to demonstrate the ability of the model to support the formulation of hypotheses and improvement of mechanistic understanding of stem cell tissue dynamics.13 Chondrogenesis of MSCs stimulated by transforming growth factor beta (TGFβ) has been modeled to determine the benefit of endogenous or exogenous TGF-β.14 Feedback regulation of tissue regeneration following injury has been mathematically described to assess how tissues regenerate following injury, while maintaining a quasi-constant fraction of stem cells.15 Stem cells can specialize into other cell types, and this phenomenon was accounted for in mathematical models by adding differentiation terms to the relevant model equations and additional equations for new compartments of cells corresponding to the modeled differentiated cell types.16 Later, the usefulness of modeling cell cycle heterogeneity was acknowledged, mainly because of its importance in research areas such as chemotherapy-based cancer treatment.17,18 Cell cycle heterogeneity can be expressed mathematically with the aid of population balance models (PBMs), which distribute the cell population over a relevant coordinate (cell property) such as size, age, or DNA content. Herein, we present a PBM that describes the osteogenic differentiation of UCB MSCs. The model utilizes heterogeneous cell cycle dynamics, osteogenic gene expression, and energy metabolism to capture stem cell differentiation toward the osteogenic lineage and to investigate the effect of changes in cell cycle duration on the differentiation efficiency of UCB MSCs. We demonstrate, by combining in silico and in vitro data, the development of distinct cell populations during osteogenic differentiation, deconvoluting the phenotypic heterogeneity of differentiating UCB MSCs.

Table 1. Model Parameters and Their Nominal Values parameter name

parameter value

dx xmin, G0G1 phase xmin, S phase xmin, G2M phase xmax, G0G1 and G2M phase xmax, S xthr σx kcat,1 kcat,2 kcat,3 kcat,4 kcat,5 kcat,6 kcat,7 kcat,8 kcat,9 kE1,1 kE1,2 kE1,3 kE2,1 kE2,2 kE2,3 kE3,1 kE3,2 kE3,3 kT1 kT2 kT3 Vcell Vreactor E(MSC) use,ATP E(Pre) use,ATP E(OB) use,ATP ωrunx ωon kTrs,runx kTrs,on kDNA,runx kDNA,on nrunx non kdecay ksRnx kson nsRnx nson

2 × 10−5 day−1 0 1 20 100 2 70 5 21.27 day−1 22.66 day−1 2.09 day−1 19.42 day−1 563.43 day−1 1081.33 day−1 1334.65 day−1 1541.34 day−1 446.09 day−1 5.25 × 10−18 L day−1 8.42 × 10−17 L day−1 1.25 × 10−18 L day−1 2 × 10−10 L day−1 2 × 10−10 L day−1 2 × 10−10 L day−1 1.47 × 10−18 L day−1 4.94 × 10−15 L day−1 1.59 × 10−17 L day−1 1000 pmol L−1 1000 pmol L−1 1000 pmol L−1 3.5 × 10−12 L 0.1 L 400 L pmol−1 80 L pmol−1 1 L pmol−1 0.2775 day−1 46.2497 day−1 4.287 day−1 600 day−1 117.0355 0.3953 15 15 1.85 day−1 0.19 35 50 50

phases G2 and M (G2M). The S (synthesis) phase was explicitly modeled. Cell cycle dynamics were described by eqs 1−10. Equation 1 was used to describe the G0G1 phase, and it contains terms for phase progress, cell death, transition to the next (synthesis) phase, transition from the previous (G2M) phase, and differentiation.

2. METHODS 2.1. Mathematical Model. The mathematical description of the osteogenic differentiation process of MSCs was split into three main partscell cycle, intracellular metabolism, and gene expressionas well as the interactions between the different modeled phenomena. Definitions of all the symbols used in the mathematical model can be found in Table 1. 2.1.1. Cell Cycle. In the description of the biological system, cells from G0 and G1 were considered aggregated into a lumped phase called G0G1 in the model equations; the same applied to

∂(kx·N (x , t )) ∂N (x , t ) + + dx ·N (x , t ) + T (x) ·N (x , t ) ∂t ∂x = δ(x − xmin) ·Ti + Di (x) − Do(x) B

(1)

DOI: 10.1021/acs.iecr.9b01988 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

where P(x) represents the normal cumulative probability function, as defined in eq 8. Note that when the differentiation fraction (s) is zero (i.e., no differentiation), eq 6 is reduced to eq 10. Cell cycle phases were described by eqs 1−10. Equation 11 links cell cycle phases via the transition rates between adjacent phases.

Cells only transition from the S phase to the G2M phase after the completion of DNA replication (at x = xmax = 2), and do not differentiate while traversing the synthesis phase. Therefore, the mathematical description of the S phase is simpler, as given by eq 2. ∂(kx·N (x , t )) ∂N (x , t ) + + dx ·N (x , t ) = δ(x − xmin)·Ti ∂t ∂x (2)

Ti(S) = To(G0 + G1)

The G2M phase was described in a similar manner to G0G1, taking into account the fact that cells do not differentiate while traversing this phase (see eq 3).

The G0G1 phase receives twice the number of cells transitioning from G2M, to account for cell division at the end of the cell cycle. 2.1.2. Metabolism. The mass balance for each of the intracellular metabolites considered by the model was expressed by eq 12. The first addend of the equation expresses how cells differentiating from type τ − 1 to type τ impact the average metabolite concentration of cells of the target type τ, while the second addend describes changes to intracellular metabolite concentrations that are due to (intracellular) metabolic reactions.

(3)

For each of the three modeled phases, the number of cells that transition to the next phase in the cell cycle was calculated using eq 4 for the G0G1 and G2M phases, and eq 5 for the S phase. In addition, each phase used its own transition function T(x) for eq 4 when computing the rate To at which cells transition to the next phase. To =

∫x

xmax

dM (j τ) dt

T (x ) · N (x , t ) d x

To = kx·N (xmax , t )

+=

∫

· N (x , t )

(7)

ri = kcat, i·Mj → ri

zz zz { dχ

rnr + i

1 − P(χ ) dχ

min

dt

(9)

ij −t ·Cells(τ) yz i zz = zz z τ = 1 k Vreactor {

∑ jjjjj 3

=

i = 1, nt

d P(x) dx

1 − P(x)

(14)

i = 1, nr

(15)

(out) kE , i M j → ti ·k T , i ti = = · Vcell Vcell k T , i + Mj → ti

dM (out) j → ti

The transition rate for the G2M phase was defined by eq 10, T (x ) = k x ·

y Np(x , t ) dx zzzz {

i = 1, nt (16)

The nutrient and waste transport across the cell membrane performed by each cell is taken into account in the culture medium, using a material balance for the extracellular metabolites:

min

∫x

xmin

The transport of metabolites across the cell membrane also contributes to the material balance of intracellular metabolites (see eq 16).

∫x 1 − P(χ ) dχ xmax

k

xmax

Equation 15 was used to describe intracellular metabolite reaction rates. The subscript notation j→ri denotes the index j of the metabolite consumed in reaction i (i.e., the reactant).

x

cpd(x) =

i

∑ jjjj∫ p=1

(8)

−∞

(13)

(6)

2y

)

Di(x) dx

min

3

Differentiation and phase transition were both shaped by the cumulative normal probability function (see eq 8). Since s was defined as the phase-distributed probability that a cell differentiates before transitioning, the cumulative probability of differentiation, cpd(x), was also defined in terms of the cumulative normal probability (see eq 9).

P(x) =

xmax

Cells(τ) =

1 − s ·cpd(x) − (1 − s) ·P(x)

(

∫x

while the total number of type τ cells (Cells(τ)) was defined by eq 14.

d

i 1 χ−x expjjjj− 2 · σ thr x k σx· 2π

STOICj , i ·ri

i=1

The term + represents the rate at which cells differentiate to type τ, and this term is defined by eq 13,

d s) dx P(x)

s · dx cpd(x)

x

∑

(12)

(5)

1 − s ·cpd(x) − (1 − s)P(x)

Do(x) = kx·

·(M (j τ) − M (j τ− 1)) + (τ ) Cells (τ )

G0G1 cells transition to the synthesis phase only if they do not first differentiate into cells of another type. For cells in G2M (unlike G0G1), exit from the phase is the only phenomenon to be described. For phase G0G1, eq 6 describes the transition rate to the synthesis phase, and eq 7 describes the differentiation rate into a more-specialized cell type: MSCs into preosteoblasts, and preosteoblasts into osteoblasts. T (x ) = k x ·

nr + nt

+(τ)

=

(4)

min

(1 −

Ti(G0 + G1) = 2·To(G2 + M) (11)

∂(kx·N (x , t )) ∂N (x , t ) + + dx ·N (x , t ) + T (x) ·N (x , t ) ∂t ∂x = δ(x − xmin) ·Ti

Ti(G2 + M) = To(S)

ji −ti·Cells(τ)

∑ jjjjj 3

τ=1

j k

Vreactor

·

y M (out) j → ti · k T , i z zz zz k T , i + Mj → ti zz {

(17)

Apart from the levels of intracellular metabolites, the model also quantifies the (intracellular) production rates of two important cofactors: ATP and NADH (see eq 18).

(10) C


Article

Industrial & Engineering Chemistry Research dEj dt

nr

=

SI(total) = Si + i

∑ STOICnm+ j ,i ·ri

j,i≠j

(18)

i=1

ne

∑ Euse,j · j=1

dEj dt

(19)

2.1.3. Gene Expression. The expression of Runx2 and osteonectin was described using eq 20 for each of the two genes and three cell types. The average (intracellular) gene expression level may change when less-differentiated cells specialize into a more-differentiated type (first addend), as it does for the intracellular metabolite concentration. The other terms in eq 20 account for the basal expression rate of the gene in MSCs, the binding of the activator to the promoter region of DNA, and the decay of mRNA, respectively. The terms +(τ) and Cells(τ) have the same meaning as for intracellular metabolism and are defined by eqs 14 and 15, respectively. Agng dg +(τ) (τ ) (τ− 1) = · − + ω + · − kdecay ·g g g k ( ) n g Trs, g ng g dt kDNA Cells(τ) g + Ag (20)

The fraction of cells that differentiate while traversing the G0G1 phase before transitioning to the synthesis phase is calculated using eqs 21. Runx2 controls the differentiation of MSCs to preosteoblasts, and Osteonectin controls the differentiation of preosteoblasts to osteoblasts. s(MSC) =

Rnx nsRnx + Rnx nsRnx

nsRnx ksRnx

s(Pre) =

(22)

The parameters were then ranked in descending order, based on their total sensitivity indices. The threshold value of 0.1 was used to determine parameter significance, regarding only the parameters with total sensitivity indices higher than the threshold for at least one of the model outputs as being significant. 2.3. Mesenchymal Stem Cell Isolation and Culture. Umbilical cord blood was used as a source of fetal mesenchymal stem cells. Mononuclear cells (MNCs) were isolated on the same day as blood collection (NHS Blood & Transplant, Colindale Blood Bank) and UCB MSCs were isolated as previously described,9,21,22 with respect to the relevant ethical approval and upon receiving informed consent from the mothers (05/Q0405/20, NRES Committee, London, Harrow). Briefly, umbilical cord blood from three donors was diluted with phosphate-buffered saline (PBS) (at a ratio of 1 volume of blood to 4 volumes of PBS) and the MNC fraction was extracted by means of Ficoll (Ficoll-paque; GE Healthcare, United Kingdom) density gradient centrifugation at 400g for 30 min without brakes, according to the manufacturer’s instructions. MNCs were washed in PBS, resuspended in basal αMEM GlutaMax-I medium (Gibco, Life Technologies, UK) with 10% fetal bovine serum (Gibco, Life Technologies, UK) and 1% penicillin/ streptomycin (Sigma−Aldrich, UK), seeded in tissue culture flasks (Corning, UK) and incubated at 37 °C, 21% O2, and 5% CO2. Two days following the initial seeding, nonadherent cells were discarded and culture of the adherent cell fraction continued until confluence. Spindle-shaped adherent MSCs were subcultured until passage 5, when osteogenic induction was initiated. Osteogenesis was performed over a 21-day period with the use of the same basal medium, supplemented with 10 mM β-glycerophosphate (Sigma−Aldrich, UK), 50 μg/mL ascorbic acid 2-phosphate (Sigma−Aldrich, UK) and dexamethasone (Sigma−Aldrich, UK) at 10−7 M, as previously described.9 2.4. Metabolomics. Metabolite concentration data were derived based on the metabolomics data set published by Klontzas et al.,9 using custom-made standard curves for all the metabolites included in the model, in a Shimadzu, Model QP2010 Ultra GC-MS system, according to previously established methods.23−25 Peak areas of glucose, pyruvate, lactate, citrate, isocitrate, succinate, fumarate, and glutamine were converted to concentration values and expressed as pM/ cell by normalizing with the concentration of ribitol (internal standard). All metabolite standards were purchased from Sigma−Aldrich UK. 2.5. qRT-PCR. RNA extraction was performed with the use of RNeasy Mini Kit, according to the manufacturer’s instructions (Qiagen, UK). Real-time quantitative RT-PCR was performed with the use of KAPA SYBR FAST One-Step ABI PRISM Kit (Sigma−Aldrich, UK), according to the manufacturer’s protocol (5 min reverse transcription at 42 °C, 3 min enzyme inactivation at 95 °C, and 40 cycles of 3 s denaturation at 95 °C and 30 s of annealing/extension at 62 °C) in a StepOne Plus qRT-PCR instrument (Applied Biosystems, CA, USA). Expression of osteonectin, and runt-related transcription factor 2 (Runx2) mRNA was quantified with the use of the following primers: osteonectin (forward) 5′-TCCACAGTACCGGATTCTCTCT-3′, (reverse) 5′-TCTATGTTAGCACCTTGTCTCCAG-3′; Runx2 (forward) 5′-GAGGTACCAGATGGGACTGTG-3′; (reverse)

The concentration of these cofactors was assumed constant, because their consumption rate was considered to be equal to their production rate. The production rate of these energy molecules then directly regulated the cell population growth rate by means of the growth terms in the population balances, as described by eq 19. kx =

∑ Si ,j

Onnson + Onnson

nson kson

(21)

2.2. Global Parameter Sensitivity Analysis. Global sensitivity analysis (GSA) was performed on the computer implementation of the mathematical model to determine which parameters influence the prediction of the model the most, given the uncertainty in the values of those parameters; the analysis in this paper used an uncertainty interval ranging from 50% to 150% of the nominal value for each parameter of the model. The sets of parameter values used for the analysis were produced using Sobol’s sequence,19 as implemented by the MATLAB function sobolset. Next, simulation data were gathered by means of repeated gPROMS simulationsone for each set of values sampled from the parameter space (Process Systems Enterprise, gPROMS, v5.1). The model outputs subjected to the analysis are the total number of cells (of all types) and the (total) number of osteoblasts, sampled at 3.5-day intervals for the first 3 weeks of the differentiation process6 values for each variable, for a total of 12 outputs. The sensitivity indices for each output variable were computed with the gui_hdmr software, which implements the random sampling−high-dimensional model representation (RS-HDMR) technique.20 For each output, the software computes a first-order sensitivity index (Si) and a secondorder contribution (Sij); from these, we computed the total sensitivity indices using eq 22. D


Article


MSCs cycled faster (higher E(MSC) use,ATP), the shift in the intracellular metabolic activity level was delayed and the effect was more pronounced when the preosteoblasts had a lower growth rate (lower E(Pre) use,ATP and higher G0G1 phase duration). Specifically, when using the nominal value of E(MSC) use,ATP, differentiation is delayed by ∼2 days (from day 8 to day 10) for E(Pre) use,ATP below the (Pre) nominal value; for high E(MSC) use,ATP and low Euse,ATP, the onset of differentiation is further delayed until approximately day 12 (see Figure 2). While all metabolite concentrations exhibited a similar fluctuation pattern within the examined time frame, the importance of glycolytic activity was different from that of the TCA cycle throughout the differentiation process, with the TCA cycle accounting for ∼87% of the total amount of energy produced by the cells during the 21-day differentiation process and glycolysis accounting for the remaining 13%. Cell cycle dynamics did not affect the onset of Runx2 expression, which always occurred between days 6 and 7 (not shown). Interestingly, an increase in growth rate (decrease in cell cycle duration) for either MSCs or preosteoblasts led to a delay in the expression of osteonectin (see Figure 2i). 3.3. Population Balance Modeling Incorporating Metabolism and Gene Switches, and Deconvoluted Heterogeneity of Differentiating Populations. The proposed mathematical model accurately captured the experimentally measured activity of glycolysis, TCA cycle, and glutaminolysis throughout the differentiation process. In-silico results indicated that levels of all described metabolites followed a similar qualitative trend from day 0 to day 21 of differentiation (Figure 3). Starting at a basal level between day 0 and day 5, the model showed a sharp concentration increase for all metabolites, which reached a maximum at day 11 before starting to steadily decrease until day 21. These results matched the experimental data, showing a global increase in the activity of glycolysis, TCA cycle, and glutaminolysis between day 6 and day 11, followed by reduced activity of the pathways at the final stages of differentiation. The mRNA levels of Runx2 and osteonectin were used to define early and late osteogenic differentiation, respectively (Figure 4), and, subsequently, the MSC, preosteoblast, and osteoblast stages. Activation patterns were successfully captured for both genes showing that Runx2 levels increased early during the differentiation period (day 3) reaching a plateau at day 13 of differentiation (Figure 4a). Osteonectin levels increased later during the differentiation period, starting approximately at day 9 and plateaued toward the final differentiation stages, at day 19 (Figure 4b). Population dynamics were captured throughout the osteogenic differentiation process, achieving deconvolution of population heterogeneity (Figure 5). Total population numbers started increasing over the first week of differentiation. The total population increased sharply between day 7 and day 14, reaching a plateau at day 14, at an approximately 4-fold higher cell number compared with the initial undifferentiated UCB MSC population. Detailed description of the individual populations appearing at every step of the differentiation process is virtually impossible by experimental means, since cells are enclosed in a dense extracellular matrix and nondestructive single cell analysis is technically complicated. Model predictions showed that the population of undifferentiated MSCs (blue line) started to decline at day 6, when osteoprogenitors first appeared. The preosteoblast (dashed orange line) population began to increase at the same time point

5′-TCGTTGAACCTTGCTACTTGG-3′; RPL13a (forward) 5′-CTATGACCAATAGGAAGAGCAACC-3′, (reverse) 5′GCAGAGTATATGACCAGGTGGAA-3′.21,26 Results were normalized with the use of RPL13a as the housekeeping gene and expressed as relative mRNA levels (2−ΔCt).27 2.6. DNA Quantification. Quantification of DNA as an indirect measure of cell numbers was performed with Quant-it PicoGreen Assay (ThermoFisher, UK) as previously published.9,28 Cells at day 7, 14, and 21 of differentiation were washed with PBS and incubated overnight in quantification buffer (1 mM EDTA, 10 mM Tris, 0.1% Triton-X (Sigma− Aldrich, UK) and 0.1 mg/mL Proteinase K (ThermoFisher, UK)). The resulting DNA solution was diluted 1:10 in PBS and 100 μL of each sample was incubated for 5 min in darkness with 100 μL of PicoGreen working solution. Fluorescence was subsequently measured in a Promega fluorescence reader (Glomax, Promega, UK) at 485 nm/530 nm excitation/ emission cell numbers were calculated using a standard curve custom-made with UCB MSCs.

3. RESULTS 3.1. Parameters Related to Gene Expression and Metabolism Affect Proliferation and Osteogenic Differentiation. Global sensitivity analysis was performed on the computer implementation of the mathematical model utilizing 32 factors (model parameters) and 2 responses (total cell count and osteoblast-only cell count). Parameters related to cell differentiation through gene expression, or to cell proliferation (Pre) through intracellular metabolism ksRnx, E(MSC) use,ATP, Euse,ATP, and kdecay, were determined to be significant (see Table 2). Gene Table 2. Significant Parameters of the Mathematical Model, As Identified Using Global Sensitivity Analysisa SI(total) i parameter name, i

for the total cell count

for the osteoblast count

kdecay kDNA,Rnx kDNA,On ω1 ksRnx ksOn kE1,1 E(MSC) use,ATP E(Pre) use,ATP kcat,2

0.119 0.069 0.137 0.142 0.272 0.093 0.115 0.194 0.162 0.104

0.215 0.105 0.103 0.087 0.109 0.108 0.060 0.069 0.059 0.066

a

Values below the threshold (0.1) are grayed out.

expression was affected by parameters such as DNA binding constants for Runx2 and osteonectin, base expression rate of Runx2 in MSCs, threshold level of osteonectin for differentiation to osteoblasts, while intracellular metabolism was influenced by parameters such as glucose uptake rate by MSCs and reaction constant for the enzymatic transformation of pyruvate into lactate (Figure 1, Table 2). 3.2. Cell Cycle Affects Metabolism and Differentiation of UCB MSCs. The model was used to assess the effect of cell cycle duration on the temporal profiles of key variables, such as the average concentrations of intracellular metabolites (for each cell type) and the levels of relative gene expression of the total cell population. Similar qualitative behavior was observed for all metabolites described by the model (Figures 2a−h). When E


Article


Figure 1. Global sensitivity analysis results. Parameters with sensitivity indices larger than the threshold at any of the time points evaluated (>0.1− dashed line) are indicated in bold red font and are deemed significant for the respective model output ((a) total cell count or (b) osteoblast-only count).

Figure 2. Influence of cell growth parameters (Euse,ATP) of mesenchymal stem cells (x-axis) and preosteoblasts (y-axis) on the differentiation onset time (z-axis), expressed as the time point at which differentiation-specific changes begin to occur for the examined variables (for average intracellular metabolite concentrations and relative gene expression levels): (a) glucose, (b) pyruvate, (c) lactate, (d) citrate, (e) isocitrate, (f) succinate, (g) fumarate, (h) glutamine, and (i) osteonectin.

4. DISCUSSION

(day 6) until approximately day 12, when the terminally differentiated osteoblasts emerged (dashed yellow line). After day 12, preosteoblast numbers decreased as a result of their terminal differentiation to osteoblasts and the number of osteoblasts increased until reaching a plateau at day 21, when their number virtually equals the total population.

Herein, the first comprehensive PBM of osteogenic differentiation has been presented, incorporating a description of cell cycle dynamics, intracellular metabolism, and gene expression. In-silico results showed that cell cycle dynamics play a significant role in the timing of the process of osteogenic differentiation of F


Article


Figure 3. Average intracellular metabolite concentrations during the osteogenic differentiation process of MSCs: (a) glucose, (b) pyruvate, (c) lactate, (d) citrate, (e) isocitrate, (f) succinate, (g) fumarate, and (h) glutamine. Comparison between the experimental values (represented as diamond symbols with vertical error bars) and the prediction of the mathematical model (solid black lines).

Figure 4. Average gene expression levels during the osteogenic differentiation process of MSCs: (a) Runx2 and (b) osteonectin. Comparison between experimental data (represented as black diamond symbols with vertical error bars) and the prediction of the mathematical model (solid black lines).

in silico during osteogenic differentiation provides a unique capability that is not achievable with conventional in vitro methods. Accounting for such heterogeneity and linking it to gene expression and metabolism can lead to culture optimization for high-quality osteogenic differentiation for bone tissue engineering purposes. The model can potentially point toward targeted optimization of culture parameters related to extracellular metabolism, leading to high-quality osteogenic differentiation. Achieving this with in vitro experimentation would require copious and costly experimentation. In order to achieve targeted culture optimization, mathematical modeling should be combined with experimentation. Such culture optimization may be achieved with the use of large-scale (e.g., genome scale) metabolic models. However, the implementation of such models is technically challenging and complicated, since stem cell differentiation demands the use of multiple objective

MSCs and confirmed that the expression of differentiation genes plays the most crucial role in the process of osteogenesis. Most importantly, the proposed mathematical model captured the heterogeneity of cell populations during osteogenic differentiation, reproducing the experimentally observed cell culture behavior over the course of osteogenesis. Stem cell culture kinetics have been previously modeled, including cell death and cell differentiation, assigning cells to compartments based on their differentiation state.29 Herein, modeling of the differentiation has been performed with higher fidelity and a greater level of detail at every step of the modeling process. Population balance enabled the quantification of individual population numbers as well as cell cycle heterogeneity. Differentiation is not only accounted for by the model, but also located precisely during the G1 phase of the division cell cycle. The ability to untangle the heterogeneity of MSC cultures G


Article


Figure 5. Partial and total cell counts during the osteogenic differentiation process of MSCs. Comparison between experimental data (black diamonds with vertical error bars) and the prediction of the mathematical model (solid black line) for the total cell count. Model predictions for the mesenchymal stem cell count (solid blue line), preosteoblast count (orange dashed line), and osteoblast count (yellow dash-dotted line) are also plotted.

Stem cells proliferate slower during differentiation,36 directing their energy through other metabolic pathways than those required for optimal (maximum) growth. In-silico experimentation affords comprehensive evaluation of cell proliferation, relative to metabolism, by linking cellular growth rate to the rate of energy production. Simulation results for parameter sweeps around their nominal value show that shorter cell-cycle durations result in delayed differentiation, even though the two cellular functions (proliferation and differentiation) do not directly compete, in terms of energy demand. This result suggests that, under the constraint of a limited time span, one may have to choose between obtaining either a higher number of less-differentiated cells or a lower number of more-differentiated cells. This effect can be explained by the fact that cell-cycle duration was strongly correlated with the duration of the G0/G1 phase; consequently, shorter cell-cycle times are translated to less time spent by cells in the G0/G1 phase and a narrower time window when differentiation can occur. Even though short cellcycle times delay the differentiation process, all the characteristic metabolic and genetic changes that we observe when using the nominal values of the parameters still occur as part of the differentiation process, but they do so at a later time. These changes include an increase in metabolic activity during the differentiation period, followed by a decrease (in metabolic activity) toward the end of differentiation, as well as an increase in the relative gene expression of Runx2, followed by osteonectin, as an indication of successful progression of osteogenesis.37 Global parameter sensitivity analysis was performed by simultaneously varying the values of all the parameters in a space ranging from 50% to 150% around their nominal values. Interestingly, analysis results show that the parameters that have the highest total sensitivity indices are directly involved in osteogenic gene expression, even though the mathematical connection between the model outputs (i.e., cell counts) and genetics first traverses the metabolic reaction rates and

functions to account for multiple competing cell behaviors throughout the differentiation process. Compared with genomescale models, the metabolic model described herein is computationally faster, since it monitors fewer metabolites; it also requires fewer measurements, only for the metabolites that were included in the mathematical model; while this mathematical model provides results for fewer metabolites (only those that it includes), it captures pathways that were shown to be significant for the osteogenic differentiation of MSCs.7,9,21,30 PBMs have been previously used to describe stem cell cultures. Wu et al. used population balance equations to describe the size of cell aggregates cultured in spinner flasks.31 Herein, population balance has been utilized in a different manner, modeling cell cycle heterogeneity as an uneven distribution of the cells between and within cell cycle phases. Our work considers cellular metabolism directly linked to cell proliferation, in comparison to Wu et al., who modeled oxygen transport as the limiting factor for stem cell proliferation. Bartolini et al. used population balance equations to describe the cell division cycle and study the proliferation of stem cells in suspension bioreactors with a special focus on bioreactor dynamics and operation, rather than intracellular gene and metabolic activity.32 A similar population balance approach has been utilized to study cell cycle phase-specific chemical treatment of leukemia,17 describing the effect of cell cycle heterogeneity on healthy and leukemic subpopulations subjected to chemotherapy. Similar to Münzer et al.,33 our work has provided a detailed description of cell cycle heterogeneity. However, our model focuses on the regulation of stem cell differentiation by gene expression and provides a more accurate description of cellular metabolism by utilizing intracellular metabolic reactions instead of extracellular metabolite measurements that may not accurately represent intracellular metabolic network activity.34,35 H



■

ACKNOWLEDGMENTS The authors of this work have received funding through the European Union Horizon 2020 research and innovation programs, under Grant Agreement No. 675585 (SymBioSys).

intracellular metabolite concentrations and is only indirectly linked with gene expression. This result is mathematically justified if we acknowledge that, in many of the simulations, the cells do not differentiate at all for certain choices for the parameter sets. In other words, while the impact of the metabolic parameters on the growth rates for individual cell types is bound at 50% around the nominal value, the values of the gene-related parameters can completely change the behavior of the model, increasing their significance value as calculated by the global analysis technique and demanding the precise determination of the values of these parameters from carefully designed experiments. This type of result requires a detailed model of gene expression during the differentiation process. Other published models of differentiation describe, in less detail, the expression of genes related to osteogenesis. While this work introduces an unstructured and segregated stem cell mathematical model to capture the cell-cycle culture heterogeneity, Vozzi et al. used a structured and unsegregated cell model in their study of the lifespan and senescent behavior of stem cells.12 Thalheim et al. employed a three-dimensional (3D) computational model to study stem cell competition in intestinal crypts involving phenomena from multiple scales, similar to this work.13 However, there are fundamental differences in the two approaches; since the intestinal crypts only contain a small number of (stem) cells (typically 5−15), they were able to model each cell individually, which would be impractical for the large cell counts that occur in osteogenic differentiation experiments (more than 50000 per well plate). The model by Chen et al. described cell numbers using two ordinary differential equations (ODEs): one for each type of cell (MSCs and chondrocytes). While their model also described MSC differentiation, they focused on the effect of the TGF-β growth factor on the chondrogenic differentiation of MSCs, while the present paper aims for the development of a mathematical model that encompasses cell cycle dynamics, intracellular metabolism, and genetics, as well as the connections between them.14 Finally, Renardy et al. investigated the stability of regenerating tissues, in terms of cell fraction by type and population recovery rate, looking for the conditions (parameter values) necessary for the tissues to maintain normal function.15 In conclusion, this work demonstrates a PBM that describes MSC osteogenic differentiation linking metabolism, gene expression, and cell-cycle dynamics. This work can enable future model-based culture optimization. Further application of the model in 3D culture systems can aid the large-scale development of clinical-grade bone grafts with orthopedic and maxillofacial applications.

■

Article

■

LIST OF ABBREVIATIONS MSCs, mesenchymal stem cells; UCB MSCs, umbilical cord blood mesenchymal stem cells; MNCs, mononuclear cells; RSHDMR, random samplinghigh-dimensional model representation; PBM, population balance model

■

REFERENCES

(1) Amini, A. R.; Laurencin, C. T.; Nukavarapu, S. P. Bone tissue engineering: recent advances and challenges. Crit Rev. Biomed Eng. 2012, 40, 363−408. (2) Giannoudis, P. V.; Dinopoulos, H.; Tsiridis, E. Bone substitutes: an update. Injury 2005, 36, S20−S27. (3) Evans, C. H. Advances in regenerative orthopedics. Mayo Clin. Proc. 2013, 88, 1323−1329. (4) Roseti, L.; Parisi, V.; Petretta, M.; Cavallo, C.; Desando, G.; Bartolotti, I.; Grigolo, B. Scaffolds for bone tissue engineering: State of the art and new perspectives. Mater. Sci. Eng., C 2017, 78, 1246−1262. (5) Klontzas, M. E.; Kenanidis, E. I.; Heliotis, M.; Tsiridis, E.; Mantalaris, A. Bone and cartilage regeneration with the use of umbilical cord mesenchymal stem cells Bone and cartilage regeneration with the use of umbilical cord mesenchymal stem cells. Expert Opin. Biol. Ther. 2015, 15, 1541−1552. (6) McLeod, C. M.; Mauck, R. L. On the origin and impact of mesenchymal stem cell heterogeneity: New insights and emerging tools for single cell analysis. Eur. Cells Mater. 2017, 34, 217−231. (7) Chen, C.-T.; Shih, Y.-R. V.; Kuo, T. K.; Lee, O. K.; Wei, Y.-H. Coordinated changes of mitochondrial biogenesis and antioxidant enzymes during osteogenic differentiation of human mesenchymal stem cells. Stem Cells 2008, 26, 960−968. (8) Shum, L. C.; White, N. S.; Mills, B. N.; de Mesy Bentley, K. L.; Eliseev, R. A. Energy Metabolism in Mesenchymal Stem Cells During Osteogenic Differentiation. Stem Cells Dev. 2016, 25, 114−122. (9) Klontzas, M. E.; Vernardis, S. I.; Heliotis, M.; Tsiridis, E.; Mantalaris, A. Metabolomics analysis of the osteogenic differentiation of umbilical cord blood mesenchymal stem cells reveals differential sensitivity to osteogenic agents. Stem Cells Dev. 2017, 26, 723−733. (10) Granchi, D.; Ochoa, G.; Leonardi, E.; Devescovi, V.; Baglìo, S. R.; Osaba, L.; Baldini, N.; Ciapetti, G. Gene expression patterns related to osteogenic differentiation of bone marrow−derived mesenchymal stem cells during ex vivo expansion. Tissue Eng., Part C 2010, 16, 511−524. (11) Drissi, H.; Hushka, D.; Aslam, F.; Nguyen, Q.; Buffone, E.; Koff, A.; Van Wijnen, A. J.; Lian, J. B.; Stein, J. L.; Stein, G. S. The cell cycle regulator p27 kip1 contributes to growth and differentiation of osteoblasts. Cancer Res. 1999, 59, 3705−3711. (12) Vozzi, G.; Lucarini, G.; Dicarlo, M.; Andreoni, C.; Salvolini, E.; Ferretti, C.; Mattioli-Belmonte, M. In vitro lifespan and senescent behaviour of human periosteal derived stem cells. Bone 2016, 88, 1−12. (13) Thalheim, T.; Buske, P.; Przybilla, J.; Rother, K.; Loeffler, M.; Galle, J. Stem cell competition in the gut: insights from multi-scale computational modelling. J. R. Soc., Interface 2016, 13, 20160218. (14) Chen, M. J.; Whiteley, J. P.; Please, C. P.; Schwab, A.; Ehlicke, F.; Waters, S. L.; Byrne, H. M. Inducing chondrogenesis in MSC/ chondrocyte co-cultures using exogenous TGF-β: a mathematical model. J. Theor. Biol. 2018, 439, 1−13. (15) Renardy, M.; Jilkine, A.; Shahriyari, L.; Chou, C. S. Control of cell fraction and population recovery during tissue regeneration in stem cell lineages. J. Theor. Biol. 2018, 445, 33−50. (16) Loeffler, M.; Wichmann, H. E. A Comprehensive Mathematical Model of Stem Cell Proliferation Which Reproduces Most of the Published Experimental Results. Cell Proliferation 1980, 13, 543−561. (17) Fuentes-Garí, M.; Misener, R.; García-Münzer, D.; Velliou, E.; Georgiadis, M. C.; Kostoglou, M.; Pistikopoulos, E. N.; Panoskaltsis,

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (A. Mantalaris). *E-mail: [email protected] (M. C. Georgiadis). ORCID

Michail E. Klontzas: 0000-0003-2731-933X Margaritis Kostoglou: 0000-0001-7955-0002 Michael C. Georgiadis: 0000-0002-2016-5131 Present Address ∇

University “Politehnica” of Bucharest, Bucharest, 060042, Romania.

Notes

The authors declare no competing financial interest. I


Article

Industrial & Engineering Chemistry Research N.; Mantalaris, A. A mathematical model of subpopulation kinetics for the deconvolution of leukaemia heterogeneity. J. R. Soc., Interface 2015, 12, 20150276. (18) Fuentes-Garí, M.; Misener, R.; Georgiadis, M. C.; Kostoglou, M.; Panoskaltsis, N.; Mantalaris, A.; Pistikopoulos, E. N. Selecting a differential equation cell cycle model for simulating leukemia treatment. Ind. Eng. Chem. Res. 2015, 54, 8847−8859. (19) Sobol’, I. M. On the distribution of points in a cube and the approximate evaluation of integrals. USSR Comp Math Math Phys. 1967, 7, 86−112. (20) Ziehn, T.; Tomlin, A. S. GUI-HDMR - A software tool for global sensitivity analysis of complex models. Env. Model. Softw. 2009, 24, 775−785. (21) Klontzas, M. E.; Reakasame, S.; Silva, R.; Morais, J. C. F.; Vernardis, S.; MacFarlane, R. J.; Heliotis, M.; Tsiridis, E.; Panoskaltsis, N.; Boccaccini, A. R.; Mantalaris, A. Oxidized alginate hydrogels with the GHK peptide enhance cord blood mesenchymal stem cell osteogenesis: A paradigm for metabolomics-based evaluation of biomaterial design. Acta Biomater. 2019, 88, 224−240. (22) Tahlawi, A.; Klontzas, M. E.; Allenby, M. C.; Morais, J. C. F.; Panoskaltsis, N.; Mantalaris, A. RGD-functionalized polyurethane scaffolds promote umbilical cord blood mesenchymal stem cell expansion and osteogenic differentiation. J. Tissue Eng. Regener. Med. 2018, DOI: 10.1002/term.2784. (23) Kanani, H.; Chrysanthopoulos, P. K.; Klapa, M. I. Standardizing GC-MS metabolomics. J. Chromatogr. B: Anal. Technol. Biomed. Life Sci. 2008, 871, 191−201. (24) Kanani, H. H.; Klapa, M. I. Data correction strategy for metabolomics analysis using gas chromatography−mass spectrometry. Metab. Eng. 2007, 9, 39−51. (25) Devito, L.; Klontzas, M. E.; Cvoro, A.; Galleu, A.; Simon, M.; Hobbs, C.; Dazzi, F.; Mantalaris, A.; Khalaf, Y.; Ilic, D. Comparison of human isogeneic Wharton’s jelly MSCs and iPSC-derived MSCs reveals differentiation-dependent metabolic responses to IFNG stimulation. Cell Death Dis. 2019, 10, 1−13. (26) Studer, D.; Lischer, S.; Jochum, W.; Ehrbar, M.; Zenobi-Wong, M.; Maniura-Weber, K. Ribosomal protein L13a as a reference gene for human bone marrow-derived mesenchymal stromal cells during expansion, adipo-, chondro-, and osteogenesis. Tissue Eng., Part C 2012, 18, 761−71. (27) Livak, K. J.; Schmittgen, T. D. Analysis of Relative Gene Expression Data Using Real-Time Quantitative PCR and the 2−ΔΔCT Method. Methods 2001, 25, 402−408. (28) Grayson, W. L.; Fröhlich, M.; Yeager, K.; Bhumiratana, S.; Chan, M. E.; Cannizzaro, C.; Wan, L. Q.; Liu, X. S.; Guo, X. E.; VunjakNovakovic, G. Engineering anatomically shaped human bone grafts. Proc. Natl. Acad. Sci. U. S. A. 2010, 107, 3299−304. (29) Deasy, B. M.; Jankowski, R. J.; Payne, T. R.; Cao, B.; Goff, J. P.; Greenberger, J. S.; Huard, J. Modeling stem cell population Growth: Incorporating terms for proliferative heterogeneity. Stem Cells 2003, 21, 536−545. (30) Pattappa, G.; Heywood, H. K.; de Bruijn, J. D.; Lee, D. A. The metabolism of human mesenchymal stem cells during proliferation and differentiation. J. Cell. Physiol. 2011, 226, 2562−2570. (31) Wu, J.; Rostami, M. R.; Cadavid Olaya, D. P.; Tzanakakis, E. S. Oxygen transport and stem cell aggregation in stirred-suspension bioreactor cultures. PLoS One 2014, 9, e102486. (32) Bartolini, E.; Manoli, H.; Costamagna, E.; Jeyaseelan, H. A.; Hamad, M.; Irhimeh, M. R.; Khademhosseini, A.; Abbas, A. Population balance modelling of stem cell culture in 3D suspension bioreactors. Chem. Eng. Res. Des. 2015, 101, 125−134. (33) García Münzer, D. G.; Kostoglou, M.; Georgiadis, M. C.; Pistikopoulos, E. N.; Mantalaris, A. Cyclin and DNA distributed cell cycle model for GS-NS0 cells. PLoS Comput. Biol. 2015, 11, e1004062. (34) Goudar, C. T.; Biener, R.; Konstantinov, K. B.; Piret, J. M. Error propagation from prime variables into specific rates and metabolic fluxes for mammalian cells in perfusion culture. Biotechnol. Prog. 2009, 25, 986−998.

(35) Villas-Bôas, S. G.; Moxley, J. F.; Akesson, M.; Stephanopoulos, G.; Nielsen, J. High-throughput metabolic state analysis: the missing link in integrated functional genomics of yeasts. Biochem. J. 2005, 388, 669−677. (36) Ruijtenberg, S.; van den Heuvel, S. Coordinating cell proliferation and differentiation: Antagonism between cell cycle regulators and cell type-specific gene expression. Cell Cycle 2016, 15, 196−212. (37) Kalajzic, I.; Staal, A.; Yang, W. P.; Wu, Y.; Johnson, S. E.; Feyen, J. H. M.; Krueger, W.; Maye, P.; Yu, F.; Zhao, Y.; Kuo, L.; Gupta, R. R.; Achenie, L. E. K.; Wang, H. W.; Shin, D. G.; Rowe, D. W. Expression profile of osteoblast lineage at defined stages of differentiation. J. Biol. Chem. 2005, 280, 24618−24626.

J


Capturing Mesenchymal Stem Cell Heterogeneity during Osteogenic

Recommend Documents