Hall et al. Journal of Biological Engineering (2017) 11:38 DOI 10.1186/s13036-017-0080-5RESEARCH Open AccessModeling the behavior of human inducedpluripotent stem cells seeded on meltelectrospun scaffoldsMeghan E. Hall1, Nima Khadem Mohtaram2, Stephanie M. Willerth2,3,4,5,6* and Roderick Edwards1,6AbstractBackground: Human induced pluripotent stem cells (hiPSCs) can form any tissue found in the body, making themattractive for regenerative medicine applications. Seeding hiPSC aggregates into biomaterial scaffolds can controltheir differentiation into specific tissue types. Here we develop and analyze a mathematical model of hiPSC aggregatebehavior when seeded on melt electrospun scaffolds with defined topography.Results: We used ordinary differential equations to model the different cellular populations (stem, progenitor,differentiated) present in our scaffolds based on experimental results and published literature. Our model successfullycaptures qualitative features of the cellular dynamics observed experimentally. We determined the optimal parametersets to maximize specific cellular populations experimentally, showing that a physiologic oxygen level (∼ 5%)increases the number of neural progenitors and differentiated neurons compared to atmospheric oxygen levels(∼ 21%) and a scaffold porosity of ∼ 63% maximizes aggregate size.Conclusions: Our mathematical model determined the key factors controlling hiPSC behavior on melt electrospunscaffolds, enabling optimization of experimental parameters.Keywords: Stem cell, Ordinary differential equation, Differentiation, Proliferation, Tissue engineeringBackgroundTissue engineering combines biomaterials and cells, cre-ating functional structures that can replace damagedregions of tissue [1, 2]. Pluripotent stem cells can differen-tiate into any cell type found in an organism, making thema valuable tool for tissue engineering [1–3]. In 2006, Taka-hashi and Yamanaka discovered that mature cells couldbe reprogrammed back into a pluripotent state, whichintroduced the option of using patient derived stem cellsas a tool for engineering tissues [4]. These cells weretermed induced pluripotent stem cells (iPSCs). Control-ling the differentiation of human iPSCs into functionaltissues remains difficult because complicated cell signal-ing networks regulate this process [5–9]. Both physicaland chemical cues control stem cell differentiation [5–9],*Correspondence: willerth@uvic.ca2Department of Mechanical Engineering, University of Victoria, Victoria, Canada3Division of Medical Sciences, University of Victoria, Victoria, CanadaFull list of author information is available at the end of the articleand the signalling mechanisms involved may be con-trolled, chaotic, random, or a combination of these types[10, 11].Stem cells can produce additional stem cells or they candifferentiate, which occurs when their gene expression isaltered to reflect a terminal state [8]. During the differ-entiation process, cells go through distinct transition to aprogenitor state where they are not considered stem cellsor differentiated cells [9, 10]. Distinct cellular markers,such as proteins, can distinguish the state of a cell in thedifferentiation process [8, 12–14]. While previously dif-ferentiation was thought to be irreversible, mature cellscan revert from the terminally differentiated state into thestem state in both nature and the lab [9, 15]. While maturecells rarely revert from a differentiated state, progenitorsrevert back into stem cells more frequently [9, 16]. Ourmodel considers these three cell states (stem, progenitor,and differentiated).These different cellular populations interact throughvarious mechanisms, including release of soluble factors© The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, andreproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to theCreative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.Hall et al. Journal of Biological Engineering (2017) 11:38 Page 2 of 20and binding of membrane proteins [17]. These mecha-nisms activate signaling pathways that vary in their speedand effective distance. Collectively these mechanisms leadto major differences in the overall effect of the signalon stem cell behavior in terms of proliferation versusdifferentiation. For example, the cellular responses to con-tact (a mechanical process) and oxygen levels (a chemicalprocess) would likely have different mechanisms. In tis-sue engineering, the scaffold provides mechanical cuesto the seeded stem cells, inducing differentiation throughmechanical input [5, 14, 18, 19]. Stiffness affects the typeof cell resulting from differentiation: cells differentiateinto neural, mesenchymal, and bone with increasing scaf-fold stiffness [19]. Moreover, substrate stiffness modulateshow neural stem cells differentiate into mature cells of thenervous system, including neurons, astrocytes, and oligo-dendrocytes [5]. In addition to stiffness, porosity affectsstem cell differentiation and proliferation because con-tact induces these processes, but excessive contact inhibitsaggregate growth. Scaffold topography also influencesstem cell growth and differentiation [14]. For example,fiber diameter of fibrous scaffolds affect aggregate growth:suboptimal fiber sizes inhibit cellular proliferation [14].Here we use melt electrospinning to fabricate our bio-material scaffolds. Melt electrospinning produces highlyreproducible engineered microfiber scaffolds with con-trollable properties, such as fiber diameter and porosity,that provide topographical cues. Our group has doneextensive work designing and fabricating such scaffoldsfor promoting the neuronal differentiation of pluripotentstem cells. While solution electrospinning is commonlyused to fabricate scaffolds for tissue engineering appli-cations, it presents many challenges such as the use oftoxic solvents and the lack of reproducibility. In either casescaffolds can degrade over time, resulting in decreasedstiffness, increased porosity, and changes in topography.Here we model how human iPSC-derived neural aggre-gates grow and differentiate on biomaterial scaffolds pro-duced by melt electrospinning poly caprolactone. Theexperimental procedure that precedes the aggregate seed-ing is outlined in Fig. 1. Our model requires the aggregatesto make contact with the electrospun substrates to initiatedifferentiation. Previous experimental observations indi-cate that neural aggregates in suspension do not progressfrom the progenitor state to the terminally differenti-ated state, but seeding them upon a biomaterial substratetriggers differentiation [14].There has been significant research into stem cell prolif-eration and differentiation, both experimentally and usingmathematical models [2, 5–7, 10, 11, 13, 20–22]. The cur-rent work incorporates multiple intrinsic and extrinsicfactors that affect stem cell population dynamics, makingit distinct from previous studies. The intrinsic character-istics include cell-cell signalling, differential responses toextrinsic effects, and state-specific metabolic properties.The extrinsic properties include scaffold effects, oxygenand waste effects, depth of culture medium, and con-trol of differentiation via growth factors. Using coarseapproximations of many of these processes allows us toinclude more experimental properties in a single model.This model gives insight into how the experimental proce-dures could be altered to maximize a specific populationof cells. In particular, we show that a physiologic oxygenlevel (∼ 5%) increases the number of neural progenitors aswell as differentiated neurons compared to atmosphericoxygen levels (∼ 21%) and a scaffold porosity of ∼ 63%maximizes aggregate size.MethodsExperimental methodsA custom-made melt electrospinning setup was usedto fabricate poly-(-caprolactone) (PCL, number averagemolar mass (Mn)∼ 45,000, Sigma Aldrich Chemical Co)biomaterial scaffolds [14]. Melt electrospinning was per-formed using nozzles with diameters of 200μm and 500 μmto fabricate scaffolds referred to as loop mesh 200 andloop mesh 500, respectively [14]. Increasing nozzle sizecorresponds to an increase in fiber diameter and adecrease in porosity [14]. The resulting porosities were23% for loop mesh 200 and 40% for the loop mesh500. Figure 2 shows two examples of the final scaffolds.hiPSCs were cultured on a Vitronectin XFTM matrixin the presence of TeSRTM-E8TM medium [23], then inSTEMdiffTM Neural Induction Medium (NIM) for 5 daysto induce neural differentiation to the progenitor cell state.Aggregates containing neural progenitor cells were thenseeded onto scaffolds and cultured in NIM for 12 daysto induce differentiation to the terminally differentiatedcell state.Bright field images were taken daily of neural aggre-gates seeded on each set of scaffolds using the IncuCyteTMlive imaging platform. The associated software was usedto measure the cell body cluster area for each of the12 day culture period. Cell viability of the neural aggre-gates seeded on these scaffolds was analyzed after 12 daysusing a LIVE/DEAD® Viability/Cytotoxicity Kit (Invitro-gen) [24, 25]. Terminal neuronal differentiation of hiPSCswas assessed by performing immunocytochemistry tar-geting the neuronal protein, Tuj1 [14, 25].Model developmentFigure 3 shows our model system and the interactionsbetween cellular populations. For the derivation of themodel, the aggregate cell population was divided intothree subpopulations: stem, progenitor and differentiatedcells, denoted respectively by S, P, and D, with the totalnumber of cells denoted by T where T = S + P + D .The feasible region for these variables is S,P,D ∈[ 0,∞).Hall et al. Journal of Biological Engineering (2017) 11:38 Page 3 of 20Fig. 1 Experimental protocol for neural aggregate formation and seeding on to melt electrospun scaffolds. Human induced pluripotent stem cells(colony shown on the left) were cultured in Aggrewell plates in the presence Neural Induction Medium (center) to form aggregates of neuralprogenitor cells. The neural aggregates are then seeded onto melt electrospun scaffolds where they differentiate into neurons (right)The rates of the cellular processes are positive (though weallow the possibility r = 0) and are given as follows:α : Death rate of stem cellsβ : Death rate of progenitor cellsγ : Death rate of differentiated cellsp1 : Proliferation rate of stem cellsp2 : Proliferation rate of progenitor cellsd1 : Differentiation rate of stem cells to progenitorcellsd2 : Differentiation rate of progenitor cells todifferentiated cellsr : Reversion rate of progenitor cells to stem cells.The units for the above rates are proportion per minute.All three subpopulations undergo death, but only thestem and progenitor populations proliferate as differenti-ation was considered terminal. Progenitor cells can revertinto stem cells.The model also includes the oxygen and waste concen-trations experienced by the aggregate, denoted by O andW, respectively, and diffusion from the air-medium inter-face. O and W are represented by percentages rangingfrom 0 to 100. The concentrations of O and W in theair surrounding the culture are denoted Oair and Wair ,with the same feasible regions as O and W. All threecell populations consume oxygen and produce waste. Theresulting model isdSdt = − αS − d1S +2p1S1 + S + P + D + rPdPdt = − βP − d2P − rP +2p2P1 + P + D + d1SdDdt = − γD + d2PdOdt = − u1S − u2P − u3D + Oflux(S,P,D,O)dWdt =w1S + w2P + w3D + Wflux(S,P,D,W ) ,(1)where the dependence of S,P, and D on O and W arethrough the rate parameters, α,β , γ , d1, d2, p1, p2, and r,and where ui are the oxygen consumption rates, wi are thewaste production rates, and Oflux and Wflux are the ratesof diffusion ofO andW between the neural aggregate andair above the medium.The units are minutes for time, centimeters for distance,and percentage for gas concentration within the medium.General structure of themodelThe “compartments” of this compartmental model consistof the populations of stem, progenitor and differentiatedFig. 2 Scanning electron microscope images of loop mesh 200 (left) and loop mesh 500 (right) scaffoldsHall et al. Journal of Biological Engineering (2017) 11:38 Page 4 of 20Fig. 3 Schematic diagram of the three cell states with cellular feedback. Black arrows indicate transitions between states. Red arrows indicatenegative feedbackcells along with the concentrations of oxygen and waste.This choice of model was based on a number of consider-ations. First, these cell states can be distinguished in thelab and cells can be held at each state. Second, each statehas unique properties, some of which have been deter-mined experimentally. Finally, the cellular scale is coarseenough that there is useful data for modeling from experi-mental work and from the literature, but fine enough thatthe results of themodel can be interpreted and transferredto the lab protocol. Each of the cell populations undergoesthe appropriate cellular processes for its state. The stemcells can proliferate, differentiate, and die. The progenitorcells undergo four processes: proliferation via division, dif-ferentiation to terminally differentiated cells, reversion tostem cells, and cell death. Differentiated cells can only die.Proliferation of earlier states is inhibited by the presenceof cells in the same and later states.Using ordinary differential equations (ODEs) to modelcell populations means that they are continuous ratherthan integer-valued variables, which can cause unrealisticresults when cell numbers fall to low values. Here the cellpopulations number in the hundreds to thousands of cells,limiting any behavioral artifacts that may arise from usingODEs, and in particular, stochastic effects.Local oxygen and waste concentrations also were incor-porated. Oxygen levels influence stem cell proliferationand differentiation [8, 20, 21, 26–28]. Additionally, O2 andCO2 influence neural stem cell differentiation and theselevels can be changed experimentally [8, 20, 21, 26–28].Thus, including these variables can help determine howto optimize current experimental protocols. The modeluses CO2 (a cellular waste product that can be measuredexperimentally) as a proxy for waste.ScaffoldmodelingThe model incorporates the scaffold properties through acell-scaffold contact rate, (C), which we take to range from1 to 10, where 1 represents a 90% porous scaffold and 10represents a solid scaffold, i.e., 0% porosity. This contactrate, C, increases with decreasing porosity, but does notgo below 1 because we consider lower values to representthe scenario where the cells would not adhere to the sub-strate. The neural aggregate is roughly spherical, so 100%contact does not mean that all the cells are in contact withthe scaffold, but rather that the maximum possible pro-portion of cells are making contact with the scaffold. Thismaximum possible contact is estimated heuristically to beabout 10%. The relationship between porosity and C isC = (100% − Porosity) · (0.1) ,where the porosity has units of percent and the factor 0.1represents the maximum possible contact.Experimentally, scaffold porosity was decreased byincreasing fiber diameter rather than by increasing thedensity of fibers with the same diameter. Altering porosityin this manner is not optimal as it also changes pore size. Itwas used in this study because of experimental limitations.We do not explicitly account for changes in fiber diam-eter and pore size in the model. However, the data usedto fit some of the effects of scaffold porosity come fromscaffolds with different fiber diameters, so the effects areimplicitly included in the model. Another factor relatedto scaffolds is the topography. The loop mesh scaffoldsare 2D scaffolds formed by randomly aligned layers ofloop fibers. We did not include the effects of topographybecause it has previously been optimized experimentally[14] and it would add a significant level of complexity tothe model.The effect of scaffold porosity on differentiation andproliferation is not linear. If the scaffold is too porous,the aggregates cannot adhere, and they fall through thescaffolds gaps. Additionally, a too-porous scaffold doesnot act on the cells strongly enough to signal for pro-liferation and differentiation. Conversely, a non-porousscaffold inhibits proliferation because contact inhibitioncomes into play. A non-porous scaffold does not affect dif-ferentiation to the same degree. The effect of contact ondifferentiation plateaus at a certain rate of contact. TheseHall et al. Journal of Biological Engineering (2017) 11:38 Page 5 of 20effects are modeled as functional terms in the prolifera-tion and differentiation rates of the stem and progenitorcell variables.Determination of functional effectsMany parameters in Eq. (1) depend on O, W, and C innon-trivial ways. These functional effects multiply thebaseline experimental rates for each of the processes.They are included multiplicatively because they are allindependent effects. Data were taken from experimentsthat were testing alteration of only one condition at atime. The qualitative and quantitative data are taken fromliterature using similar cells and under similar cultureconditions in order to minimize differences arising fromfactors that were not of interest.Each effect was first determined qualitatively, then fitto quantitative data. The functions were fitted manuallybecause limited data were available. All of the functionsfor the effects take values greater than 0, where valuesabove 1 increase the rate and values below 1 decrease therate from the baseline value. Once the functional effectsare determined, it is useful for later analysis to find themaximum andminimum values for each. Table 1 gives theindividual functional effects, as well as the ranges for eachover the appropriate domains of O, W and C. Figure 4shows the individual effects and Table 2 shows how thesefunctions are built into the parameters.The complete parameter coefficients of the model areproducts of the above functional components and theexperimental parameter values.Determination of experimental ratesThe experimental rates were determined from multi-ple sources, including experimental data and data fromliterature based on similarity to our experimental set-up.Considerations include cell type and culture conditions(O2 concentration of 21%, CO2 concentration of 5%, tem-perature of 37 ◦C, etc.).For each parameter, a combination of the functions inTable 1 multiplies the experimental parameter value toproduce the compound parameter value as detailed inTable 2.In terms of notation, the compound parameter is givenby α, and the experimental parameter value is referred toas α¯, and similarly for the other parameters.DiffusionDiffusion through the medium from the air to the neu-ral aggregate, and from the neural aggregate to the air,leads to changes in the concentrations of oxygen andwastearound the neural aggregate during the experiment. Thediffusion of oxygen and waste is affected by the scaffoldporosity. Lower porosity scaffolds decrease local diffu-sion by limiting flow under the aggregate. This effect ismodeled by terms multiplying the diffusion terms, butthe form of this functional effect could only be estimated,because quantitative data were not available.Oxygen and waste equations The complete oxygen andwaste equations are comprised of consumption and pro-duction of oxygen and waste, respectively, by each cellpopulation, and diffusion, occurring in or out depend-ing on the relative concentrations of gases between thelocal environment, i.e.,O andW, and the external environ-ment, i.e., Oair and Wair . The rates for consumption aretaken from literature [29] and converted to a percentageas previously described. Each cell population has a specificconsumption rate because metabolic requirements varywith differentiation state. Stem cells consume high levelsof oxygen [29, 30] as they complete the cell cycle fasterTable 1 Components of functional effects on parameters and extremal valuesFunction Minimum point Maximum point Function rangef1 11+10O − 0.3Oe−0.3O + 1 (3.57, 0.66) (0, 2) [ 0.66, 2]f2 11+5O − 0.3Oe−0.3O + 1 (3.77, 0.69) (0, 2) [ 0.69, 2]f3 11+2O − 0.5Oe−0.3O + 1 (3.89, 0.51) (0, 2) [ 0.51, 2]f4 (0.4O − 1)(0.4O + 1)e−0.5O + 1.01 (0, 0.01) (5.2, 1.26) [ 0.01, 1.257]f5 (0.6O − 1)(0.6O + 1)e−0.5O + 1.01 (0, 0.01) (4.6, 1.67) [ 0.01, 1.67]f6 109(1 − 11+e−0.1(W−25))(100, 0) (0, 1.03) [ 0, 1.03]f7 1 + 11+e−0.5(W−16) (0, 1) (100, 2) [ 1, 2]f8 15O+5 (100, 0.14) (0, 3) [ 0.14, 3]f9 CC+1 (1, 0.5) (100, 0.99) [ 0.5, 0.99]f10 4C3e−0.5−0.8C (100, 0) (3.75, 6.37) [ 0, 6.37]f11 200200+C (10, 0.95) (1, 0.99) [ 0.95, 0.99]f12 100100+C (10, 0.91) (1, 0.99) [ 0.91, 0.99]f13 11+e−10(O−0.25) (0, 0.076) (100, 1.00) [ 0.076, 1.00]Hall et al. Journal of Biological Engineering (2017) 11:38 Page 6 of 200 5 10 15 20 25 30 35 40O00.511.522.53Function Valuef1f2f3f4f5f80 10 20 30 40 50 60 70 80 90 100W00.511.522.5Function Valuef6f71 2 3 4 5 6 7 8 9 10C01234567Function Valuef9f10f11f12Fig. 4 Functional Effects of O (top),W (middle), and C (bottom). When applicable, black markers indicate experimental data used for fittingthan their differentiated counterparts [7, 31, 32]. As cellsdifferentiate, the cell cycle lengthens, and slower cyclingtimes lead to decreased oxygen consumption [7, 31, 32].Note that the oxygen consumption rates are multiplied byf13 (see Table 1) tomodel how cells alter oxygen use duringanoxia.Cells also produce waste products, which affect cellu-lar processes upon accumulation [33, 34]. Thus, cells mayslow or arrest proliferation, and even trigger apoptosis ina high waste environment [33, 34]. In this model, we useCO2 as a proxy for cellular waste. According to literature,the production of CO2 closely matches the consumptionof O2 (units in mol). Thus, the values for O2 consumptionwere used to determine the CO2 production, using adifferent factor for the mol to % conversion, which wasdetermined in a similar fashion to the oxygen conver-sion factor. The waste production rates are multiplied byf13 (see Table 1) to model the way in which cells altermetabolic activity during anoxia.Diffusion equations Typical cell culture medium has arelative density to water of 1.00−1.06 (STEMCELL Tech-nologies, personal communication, June 9, 2016). Thus,we used the density of water (1.00) for our model. Thediffusion terms are given byOflux =JO2(S,P,D,O) · Aagg(S,P,D) ,Wflux =JCO2(S,P,D,W ) · Aagg(S,P,D) ,Hall et al. Journal of Biological Engineering (2017) 11:38 Page 7 of 20Table 2 Experimental and compound parameter valuesPar Experimental value range (×10−5) Ref Functional effect range Compound value rangeα = α¯f1f7 [ 0.1, 2.6]∗ [14] [ 0.66, 4] [ 0.00000066, 0.00010]β = β¯f2f7 [ 1.6, 2.6] [14] [ 0.685, 4] [ 0.000011, 0.000104]γ = γ¯ f3f7 [ 1.6, 2.6] [14] [ 0.51, 4] [ 0.00000816, 0.000104]p1 = p¯1f5f6f10 [ 69, 120] [7, 29, 37] [ 0, 11.81] [ 0, 0.014]p2 = p¯2f5f6f10 [ 45, 160] [7, 29, 37] [ 0, 11.81] [ 0, 0.019]d1 = d¯1f4f9 [ 10, 17] [8] [ 0.0050, 1.24] [ 0.00000050, 0.00021]d2 = d¯2f4f9 [ 7.3, 8.2] [8] [ 0.0050, 1.24] [ 0.00000036, 0.000102]r = r¯f8 [ 0.1, 17]∗ [8] [ 0.14, 3] [ 0.00000014, 0.00051]The functional effect is the product of feedbacks for each parameter, e.g. fα = f1f7. The compound value is the experimental value, denoted by a bar above the parameter,multiplied by the functional effect, e.g. α¯fα*No measurements available. Closest related measurements were taken. For α, the upper bound for β was used. For r, the upper bound for d1 was used. In both cases, it istaken that 0.000001 is the lower boundwhere J is the flux and Aagg is the cross-sectional areaof the aggregate (see equations in the Appendix). Notethat JO2 and JCO2 are also dependent on the volume ofmedia used and external concentrations of O2 and CO2,respectively, which are constant for any given experimen-tal setup. The flux terms forO andW are multiplied by f10and f11 (see Table 1), respectively, to include the effect ofscaffold density on diffusion.Fixed point existence and stabilityAnalyzing the dynamics of the model gives useful infor-mation that can be applied to experimental work. First, wewill show that the model predicts convergence to a steadystate (fixed point), the final resting state of the growthprocess.A fixed point of the model corresponds to a special statein the cell culture that once reached would be maintained,an aggregate with constant size and constant populationsof the three cell types. A main goal is to obtain an aggre-gate with as large a population as possible of a particularcell state. Thus, we wish to maximize the fixed points withrespect to a certain cell state. Although ultimately the goalis to produce as large a population, D, of differentiatedcells as possible, one approach might be to optimize S orP first, and then introduce a chemical factor to trigger dif-ferentiation. In order to achieve an optimal population ofone of the cell states, we can alterO,W, C, and differentia-tion rates, which are accessible by changes to experimentalprocedures. O and W can be controlled by altering theO2 and CO2 levels in the culture chamber. The value ofC can be altered by changing the porosity of the scaffoldon which the cells are seeded. The differentiation ratescan be modified by the addition of chemical factors in themedium, either at the start of the experiment, or duringthe culture period. In the following sections we considerseparately optimization of each population and total cellnumbers.Possible changes to experimental procedureIn addition to scaffold porosity (contact parameter) andseeding protocol (initial populations of each cell type),oxygen serves as a critical factor for experiments affect-ing stem cell culture. The current experimental procedureuses an oxygen controlled chamber that holds the con-centration near ambient levels (∼ 21%). Changing theoxygen concentration experienced by the cells affects bothdifferentiation and proliferation [8, 27]. In the case of neu-ral stem cells, the oxygen concentrations in the brain arelower than in the rest of the body, reaching as low as0.55% in some regions of the brain. Other areas of thebody maintain up to 9% oxygen [20]. Changing the cellculture medium influences cell behavior. Flow chamberscan continuously refresh the medium with oxygen andnutrients and remove excess waste products. We deter-mined whether a static culture or a culture with a mediumreplacement regimen would be optimal by running themodel under the conditions where oxygen and waste aretaken as variables or held constant.Relation of model dimension to experimental conditionsThe main goal of this work is to construct and analyze a5-dimensional model, but analysis of lower dimensionalversions of the model allows for a comparison of differ-ent experimental procedures. In our experimental work,we attempted to seed with an aggregate consisting onlyof progenitor cells (so that initially S(0) = D(0) = 0).If reversion of progenitors to the stem cell state is negli-gible (r = 0), then the stem cell population remains atzero, and can be excluded from the model. However, weare interested in the effects of reversion and the possibil-ity that some of the initial cells are stem cells. Oxygen andwaste can be kept at a constant level experimentally bycontinuous medium replacement. We therefore analyzedour system under four conditions, each with a different set(and number) of variables:Hall et al. Journal of Biological Engineering (2017) 11:38 Page 8 of 20• The PD system refers to the system without stemcells and with O and W kept constant by mediumreplacement;• The PDOW system refers to the system withoutstem cells and with variable O and W (no mediumreplacement);• The SPD system refers to the system with stem cellsand with O and W kept constant by mediumreplacement;• The SPDOW system refers to the system with stemcells and with variable O and W (no mediumreplacement).For the sake of brevity, the systems without stem cells arediscussed in the Appendix.ResultsExperimental resultsNeural aggregates derived from hiPSCs were seeded ontwo different scaffolds for 12 days. The average results of3 experiments 12 days after seeding are summarized inTable 3.Figure 5 and Table 4 show that both scaffold topogra-phies support cell adhesion and cell migration.As shown in Tables 3 and 4, cell body cluster area ofneural aggregates cultured on more porous scaffolds wasconsistently larger than that of neural progenitors seededon less porous scaffolds.Model analysis resultsThe SPD system is comprised of the first three equationsof System 1, without the differential equations for O andW. The corresponding schematic is given in Fig. 3.Our strategy is as follows. First, we find the fixed pointsof the system in terms of the parameters, which we doby means of an approximation, and then determine theirstability to establish conditions under which a positivefixed point is reached. Since the fixed points depend onthe parameters, we can determine how each population isaffected by each of the parameters individually, and findthe parameter values for which the population is maxi-mized. However, the parameters are not independent inexperiments, but depend onO,W, and C, so we next opti-mize the parameters as functions of O, W, and C, guidedby the independently optimized parameter sets. Using theoptimal O and W, we determine the corresponding OairTable 3 Comparison of data from three experiments for twoscaffold porosities 12 days after seedingScaffold type Loop mesh 200 Loop mesh 500Porosity (%) 40 23Tuj1 fluorescence (%) 71.5 ± 1 58.4 ± 3Cell body cluster area (mm2) 2.04 ± 0.1 0.87 ± 0.27andWair in the SPDOW model, indicating optimal cultureconditions.Existence of fixed pointsThe fixed points of the 3D system are determined byfinding the intersections of the nullclines, given by0 = − αS + rP + 2p1S1 + S + P + D − d1S0 = − βP − rP + 2p2P1 + P + D + d1S − d2P0 = − γD + d2P .(S,P,D) = (0, 0, 0) is one solution to this system. Tofind a non-zero solution, the three nullclines are solvedsimultaneously. From the D nullcline we haveD =d2Pγ,which can be substituted into the equations for the S and Pnullclines. As D is only dependent on P, finding the inter-section of the nullclines for S and P will be sufficient tofind the fixed point. The P nullcline, denotedN1, becomesP2(β + r + d2d1)− SP − S(γγ + d2)+ P(γ (β + r + d2 − 2p2)d1(γ + d2))= 0 ,or, equivalently,S = P(a − bc + P),where a = β+r+d2d1 , c =γγ+d2 and b =2p2cd1 .Similarly, the S nullcline, denoted N2, becomesS2(α + d1) + SP((α + d1)(γ + d2)γ− r)−P2( r(γ + d2)γ)+ S(α + d1 − 2p1) − rP = 0 .Therefore, the fixed points of the system are given by theintersections of two hyperbolas, and there can be betweenzero and four in principle. As (S,P,D) = (0, 0, 0) is aknown solution, there is at least one intersection, and atmost three positive intersections. It can be shown thatthe nullclines here have a unique positive intersection, theother two being negative in either S or P (and therefore notmeaningful). The S value of this positive intersection is anincreasing function of the P values. Solving the nullclinesfor a positive fixed point directly becomes too compli-cated to be useful analytically, so we approximate usingthe asymptotes of the nullclines for large values of P.The asymptotes of N1 and N2, denoted S1A and S2Arespectively, are determined by ensuring that limP→∞ S−Hall et al. Journal of Biological Engineering (2017) 11:38 Page 9 of 20a bc d ef g hFig. 5 Experimental results. Top: Fluorescence images of neuronal marker Tuj1 expressed in neural aggregates 12 days after seeding on loop mesh200 (a) and loop mesh 500 (b) scaffolds. Scale bar is 400 μm.Middle: Bright field images of neural aggregates on loop mesh 200 scaffolds 0 (c), 6 (d)and 12 (e) days after seeding . Bottom: Bright field images of neural aggregates on loop mesh 500 scaffolds 0 (f), 6 (g) and 12 (h) days after seedingS1A = 0, where S is the value of S on N1, and similarly,limP→∞ S − S2A = 0, where S is the value of S on N2. Theresulting asymptotes areS1A =P(β + r + d2d1)− 2p2γd1(d2 + γ ) ,S2A = Prα + d1 +2p1rγ(α + d1)[ (α + d1)(γ + d2) + rγ ] .Table 4 Cell body cluster area for two neural aggregates seededon loop mesh 200 and loop mesh 500Loop mesh 200 Loop mesh 500Cell body Number Cell body NumberDay cluster area (mm2) of cells cluster area (mm2) of cells0 0.85 5066 0.75 41992 0.92 5706 0.78 44544 1.11 7561 0.80 46266 1.44 11,172 0.88 53378 1.67 13,953 0.91 561210 1.82 15,874 1.00 645612 2.24 21,675 1.1 7459Number of cells was calculated from cell body cluster area assuming an initialpopulation of 4500 cells (see Appendix)Figure 6 shows an example of the fixed point givenby nullcline intersection and the approximate fixed pointdetermined by the asymptote intersection.The S, P, D, and T values at the positive asymptoteintersection areS∗ =2γ r[p1(β + d2 + r) + p2(α + d1 + rγγ+d2)][ (α + d1)(β + d2) + rα] [ (α + d1)(γ + d2) + rγ ] ,(2)P∗ =2γ[p1rd1 + p2(α + d1)2 + p2rγ (α+d1)γ+d2][ (α + d1)(β + d2) + rα] [ (α + d1)(γ + d2) + rγ ] ,(3)D∗ =2d2[p1rd1 + p2(α + d1)2 + p2rγ (α+d1)γ+d2][ (α + d1)(β + d2) + rα] [ (α + d1)(γ + d2) + rγ ] ,(4)T∗ = 2p1[rγ (β + d2) + rd1(γ + d2) + r2γ][ (α + d1)(β + d2) + rα] [ (α + d1)(γ + d2) + rγ ]+2p2[(α + d1)((α + d1)(γ + d2) + 2rγ ) + r2γ 2γ+d2][ (α + d1)(β + d2) + rα] [ (α + d1)(γ + d2) + rγ ] .(5)Our fixed points should occur at relatively large val-ues of P and S, making these approximate values accu-rate. Since the approximate fixed point values have beenHall et al. Journal of Biological Engineering (2017) 11:38 Page 10 of 20-50 0 50 100 150 200P-1000010002000300040005000S N1N2Asymptote of N1Asymptote of N2Fig. 6 Example of hyperbola and asymptote intersection for N1 and N2. Note that the intersection of the asymptotes is near the intersection of thehyperbolas at large Pexpressed explicitly in terms of the parameters, the effectsof parameters on the fixed points can be determined, andeach population can be optimized with respect to eachof the parameters, as was done for the PD system. Again,optimization of the individual populations as well as thetotal population may be of interest, depending upon thegoal of the experimentalist.Stability of fixed pointsThe Jacobian of the SPD system is given byJ(S,P,D) =⎡⎢⎣2p1(1+P+D)(1+S+P+D)2 − (α + d1) r − 2p1S(1+S+P+D)2 −2p1S(1+S+P+D)2d1 2p2(1+D)(1+P+D)2 − (β + d2 + r) −2p2P(1+P+D)20 d2 −γ⎤⎥⎦ .By the Routh-Hurwitz condition, the zero fixed point isstable when1. 2p1 + 2p2 < α + β + d1 + d2 + r and2. (2p1 − α − d1)(2p2 − β − d2) > (2p1 − α)r.and unstable when either of the inequalities is reversed.This implies that the zero fixed point is stable when either1. 2p1 < α and 2p2 < β + d2, or2. 2p1 < α and 2p2 > β + d2 andr > (2p2 − β − d2)(d1α−2p1 + 1), or3. α < 2p1 < α + d1 and 2p2 < β + d2 andr < (β + d2 − 2p2)(d12p1−α − 1),and not otherwise (neglecting equalities). Thus, the zerofixed point can only be stable if the two proliferation ratesare low. It is always stable if the proliferation rate of stemcells is small in relation to their death rate and the pro-liferation rate of progenitor cells is small in relation tothe combined rate of loss of progenitors to death and dif-ferentiation. It can still be stable when one of the twoproliferation rates is higher, but only under other restric-tions. If the proliferation rate of progenitors is higher, thenthe reversion rate of progenitors to stem cells must alsobe sufficiently large. On the other hand, if the prolifera-tion rate of stem cells is higher than half their death ratebut not as high as half the combined rate of loss of stemcells to death and differentiation, then the zero fixed pointis stable only if the reversion rate is sufficiently small. Allcases are possible within the parameter ranges given inTable 2. Thus, a stable zero fixed point is possible withinour parameter space, though it is unstable at optimalparameter values found below. The stability of the positiveequilibrium was determined numerically at the optimalparameter values found below. In all cases it proves to bestable.Optimizing the positive SPD fixed pointThe effects of parameters on the values of the positivefixed point are analyzed. We identified parameter valuesfor maximizing each subpopulation.Maximizing S∗Equation (2) shows immediately that S∗ is maximized atmaximum allowed values of p1 and p2. Differentiating S∗with respect to each of the other parameters shows thatthe maximal S∗ is achieved when α, β , d1, and d2 areminimized, and when γ and r are maximized.Maximizing P∗ and D∗Equation (3) shows immediately that P∗ is maximizedwhen p1 and p2 are maximized, and when β and d2 areminimized. Differentiating P∗ with respect to each of theother parameters shows that the maximal P∗ is achievedwhen α and r are minimized (but see below for r), when γis maximized, and when d1 = d∗1, an intermediate value inthe allowed range for d1. When the other parameters areat their optimal values, d∗1 = 0.0000008675.Hall et al. Journal of Biological Engineering (2017) 11:38 Page 11 of 20The only parameters that differ inD∗ from P∗ are d2 andγ , so the rest of the optimal parameters remain the samefor D∗ as for P∗. Differentiation with respect to d2 andγ shows that P∗ is maximized when γ is minimized, andwhen d2 = d∗2, an intermediate value in the allowed rangefor d2. When the other parameters are at their optimalvalues, d∗2 = 0.000009495.Taking the minimum allowed r maximizes P∗ and D∗.This conclusion depends on the range for r, [ rmin, rmax].Considering P∗ (and thus D∗) as a function of r, thereis a maximum at a positive r value. If we use the previ-ously determined optimal parameter values and considerour ranges for d1 and d2, then r∗ < rmin where rmin istaken from Table 1. If we had rmin < r∗, then r∗ would bethe optimal value rather than rmin. This is an importantconsideration given that the parameter ranges were deter-mined using limited data andmay not fully reflect the truevalues.Maximizing T∗From Eq. (5), the optimal p1 and p2 for maximizing T∗are their maximum values. The sign of the derivative ofT∗ with respect to each other parameter determines that,within the given parameter ranges, T∗ is maximized whenα, β , d1, and d2 are minimized, and when γ and r aremaximized. The conclusions for d1 and r are valid at leastwhen β > α and γ > α, which is the case when all otherparameters are at their optimal values.It is notable that the optimal parameter values to maxi-mize the total cell population at equilibrium are the sameas those that maximize the stem cell population. However,there is an important difference in the conditions for thisresult. This parameter set is the optimal set for all positiveparameter values when maximizing S∗, while the optimalvalues of r and d1 are dependent on the relative values ofα, β and γ when maximizing T∗. If the death rate of stemcells, α, exceeds the death rates of progenitor and differ-entiated cells, β and γ respectively, the optimal values ofr and d1 are not the same as for maximizing S∗. There-fore, the total cell population and the stem cell populationare optimized together, unless the death rate of stem cellssurpasses the death rates of the other compartments.Summary of optimization resultsFor each subpopulation, the optimal parameters areshown in Table 5. The entries with max and min indi-cate that, for that population, the optimal value for thatparameter is the maximum or minimum.Optimization with respect to oxygen, waste and contactThe above fixed points are optimized assuming thatthe parameters can be adjusted independently of eachother. This cannot be achieved in experiment because theparameters actually depend on O, W and C. In order toTable 5 Summary of optimization results for the SPDmodelParameter Max S∗ Max P∗ Max D∗ Max T∗α Min Min Min Minβ Min Min Min Minγ Max Max Min Maxp1 Max Max Max Maxp2 Max Max Max Maxd1 Min d∗a1 d∗a1 Minad2 Min Min d∗a2 Minr Max Mina Min MaxaNote that the parameter sets for max S∗ and max T∗ are the sameaIndicates a condition on the optimal valueprovide useful feedback for experimental protocols, thisdependence must be taken into account. Thus, the func-tional effects of O,W and C on the parameters were opti-mized by determining the values of O,W and C that mostclosely result in the optimal parameter set. For example,the value of p2 should be maximized in all cases, so theoptimal value of each functional effect on p2, i.e., f5, f6,and f10, should be maximized. To maximize γ , one has tomaximize f3 and f7, but there is a trade-off in attempting tomaximize f6 and f7 with respect toW simultaneously. Val-ues of W and C were estimated from the function graphsto optimally resolve these trade-offs where they occur. Forthe optimization of each of P∗,D∗, and T∗, the values usedwere W = 5 and C = 3.75. In principle, W = 0 is opti-mal, but W = 5 was taken because this level of CO2 alsoacts as a buffer in the medium [35] and does not greatlyaffect the functional values compared toW = 0. Note thatC = 3.75 agrees in a general sense with the experimentalresults on porosity of the scaffold.Optimal O values cannot easily be read from the func-tion graphs, and had to be determined by numericallyfinding the value of O for which each of the equilib-rium values, S∗, P∗, D∗ and T∗, takes its maximal value,using the already-determined optimal values of C and W.The resulting O is the optimal oxygen level that shouldbe used in the experimental protocol for maximizing theappropriate equilibrium population.The relevant functions for optimizing β and p2, namelyf2 and f5, as well as the relevant function for minimiz-ing γ , in order to maximize D∗, namely f3, have optimalvalues of O grouped near a common oxygen value. Therelevant function for optimizing d2, namely f4, and the rel-evant function for maximizing γ , in order to maximize P∗,have zero as the optimal value ofO. As d2 takes its optimalvalue at a different oxygen level than other parameters,it is of interest to consider d2 separately when determin-ing the optimal O as discussed above. Thus, d2 was takenas a free parameter in the model, with the other parame-ters remaining functions of O. Treating d2 independentlyHall et al. Journal of Biological Engineering (2017) 11:38 Page 12 of 20not only increases the equilibrium population of interest,it indicates how critical the value of d2 is to the dynam-ics of the system. In terms of experimental procedure, thisseparation of d2 from O is feasible. The differentiationprocess represented by d2 can be modified chemically.Thus, the results from this decoupling of d2 from O couldbe transferred to the experimental procedure.For similar reasons, we also considered the effect ofallowing d1 and r to be independent, to reflect experimen-tal control of the differentiation rate of stem cells and thereversion rate of progenitors to stem cells.The fixed point values for the optimized populations aregiven in Tables 6. The parameters used in calculating thesevalues are the experimental and compound values as givenin Table 2. For the dependent parameters, the appropriateoptimalO,W, and C for each situation were applied to theexperimental parameter values to give the final compoundparameter values.Including variable oxygen andwaste in themodel with stemcellsThe SPDOW model is given in System (1). The additionof the O and W variables makes a full analysis impos-sible without using numerical simulations. However, wecan use the results obtained in the SPD system to drawconclusions about the SPDOW system.In the SPD system, we determined the optimal popula-tions and the associated parameters. Since those param-eters are functions of oxygen and waste, optimal oxygenand waste concentrations can be calculated from the opti-mal parameter values. Using the approximate SPD fixedpoint values for S, P, and D, as well as the optimal O andW values, the SPDOW fixed point equations for O andWcan be numerically solved for Oair and Wair . It can easilybe shown thatO andW must then converge to these fixedpoint values. Then, with the resulting set of parametervalues, the true fixed points, rather than the approxi-mate fixed points determined by the nullcline asymptoteintersection, of the SPDOW system can be determinednumerically. The resulting SPDOW fixed points are givenin Table 6. The Oair and Wair values required to achievethe optimal O andW at the equilibrium values of S, P andD are also given in Table 6. The differences between pop-ulation sizes when various sets of parameters are madeindependent is further illustrated by the simulations inFig. 7.Effect of porosity on growth and differentiationScaffold porosity affects not only differentiation, but alsothe growth of neural aggregates seeded on the scaffold.As shown in Tables 3 and 4, there is a significant increasein the neural aggregate size when the scaffold porosity isincreased from 23 to 40%. To compare the experimen-tal results to the model results, we ran simulations of themodel withC = 7.7 (23% porosity) andC = 6 (40% poros-ity). Results are shown in Fig. 8. The dependent parameterset and the experimental parameter set for optimal D∗were used.Oxygen levelsIt has also been observed experimentally that culturingneural stem cells in physiologic oxygen conditions (∼ 5%)rather than ambient oxygen conditions (∼ 21%) increasesTable 6 Summary of SPDOW optimized populationsOptimized population S∗ P∗ D∗ T∗ O∗ W∗ Oair WairS∗ , indep 96869.0 165.2 0.6 97034.8 NA NA NA NAS¯∗ , dep 3383.3 179.2 94.3 3656.8 0.5994 5 1.9622 3.3338Sˆ∗ , dep 3385.2 179.2 94.3 3658.7 0.5993 5.0002 1.9622 3.3338S¯∗ , d1, d2 indep 81249.0 252.3 4.7 81506.0 4.1506 5 6.1227 2.5887Sˆ∗ , d1, d2 indep 81110.0 251.6 4.7 81366.3 4.1499 5.0008 6.1227 2.5887S¯∗ , d1, d2, r indep 83213.0 142.3 2.7 83358.0 4.2405 5 6.1682 2.5990Sˆ∗ , d1, d2, r indep 82908.0 141.6 2.6 83052.2 4.2029 5.0020 6.1682 2.5990P∗ , indep 1896.3 3428.2 12.0 5336.5 NA NA NA NAP¯∗ , dep 12.8 147.3 120.4 280.5 0.8544 5 1.2588 4.5059Pˆ∗ , dep 19.2 153.3 123.0 295.5 0.8419 5.0152 1.2588 4.5059P¯∗ , d2 indep 37.4 2880.5 53.7 2971.6 4.2776 5 4.9933 4.1258Pˆ∗ , d2 indep 37.5 2881.2 53.7 2972.4 4.2776 5.0001 4.9933 4.1258D∗ , indep 848.4 883.9 1028.5 2760.8 NA NA NA NAD¯∗ , dep 1.0 54.1 334.7 389.8 3.9716 5 4.4258 4.4450Dˆ∗ , dep 1.3 54.4 336.7 392.4 3.9704 5.0014 4.4258 4.4450Numerically calculated values of true equilibria are denoted Sˆ∗ , Pˆ∗ , and Dˆ∗ , while approximate values based on asymptotes are denoted S¯∗ , P¯∗ , and D¯∗Note that the independent P∗ and D∗ equilibria are unstable; all other equilibria are stableHall et al. Journal of Biological Engineering (2017) 11:38 Page 13 of 200 1 2 3 4 5 6 7 8 9 10Time (min) 10610-2100102104106Cell NumberS P D TFig. 7 Example of effects of dependence versus independence of parameters on populations. All parameters independent (dot); d1, d2, r independent(solid); d1, d2 independent (dash); all parameters dependent (dash-dot). Parameter set used is for maximizing S∗ , with C = 3.75 and O∗ ,W∗ , Oair , andWair values from Table 6. Experimental parameters (×10−5): α¯ = 0.1, β¯ = 1.6, γ¯ = 2.6, p¯1 = 119, p¯2 = 160, d¯1 = 10, d¯2 = 7.29, r¯ = 17. Independentparameters (×10−5): d1 = 0.05, d2 = 0.03645, r = 51both cell proliferation and level of differentiation [8].Because these experiments often used solid substrates forseeding, simulations were run with C = 10 (0% porosity)with dependent parameters and the experimental param-eter set for maximalD∗ to determine if the model capturesthese experimental results (Fig. 8). The simulations showthat both differentiation and growth increase when oxy-gen is at the physiologic level of 5% versus the ambientoxygen level of 21%.Timing of growth and differentiationThe timing of differentiation and growth was also of inter-est. To determine whether it would be better to grow alarge stem cell population, then instigate differentiation,rather than pushing differentiation at seeding, we ran asimulation starting at the S∗-optimized equilibrium withparameters set to optimize D∗. The results are shown inFig. 8. D rises to a much higher level than in the straight-forward D-optimizing strategy, but then falls again tosettle at an equilibrium value less than this peak level.DiscussionIn this work, we have developed a mathematical model ofhuman iPSC proliferation and differentiation on scaffolds.In an effort to derive a realistic model, the model devel-opment included both qualitative and quantitative datafrom both laboratory research and previous literature.We included variables and parameters that would notonly make interpretation of mathematical and numericalresults for experimental procedure possible, but also inorder to make optimization of desired outcomes math-ematically feasible. By analyzing the model, we are ableto explore the effects of changes in parameters thatdepend on experimental protocol, and thus guide futureexperiments. In some cases the effects are not intuitivelyobvious, as there are trade-offs between positive and neg-ative effects of change in even a single parameter.Computer-aided analysis showed that the model hasa unique positive steady state, (S∗,P∗,D∗), to which thesystem converges over time from any initial set of pop-ulation values (apart from the inevitable steady state atzero — no cells present). The stem, progenitor, and dif-ferentiated cell populations at steady state were optimizedwith respect to each individual parameter with resultsgiven in Table 5. Many of the optimization results agreewith intuitive expectations: it makes sense to maximizeproliferation and minimize the death rate of stem or pro-genitor cells in order to maximize their respective equi-librium population sizes. The model analysis also revealedthat certain parameters, such as differentiation rates, donot always have such intuitively obvious optimal values.The dependence of the model parameters on oxygen,waste and contact rate imposes additional practical con-straints on the accessibility of optimal parameter values inactual experiments. Including these constraints and waysof loosening some constraints by experimental procedure(such as introduction of differentiation factors into theculture medium) allows the mathematical model to give arealistic optimal procedure.Maximizing S∗One possible objective of an experiment is to maximizethe stem cell population, possibly with the idea of laterintroducing a differentiation factor to drive differentia-tion. In this case, the model confirms our expectation thatthe equilibrium stem cell population is maximized at theminimum value of the stem cell death rate (α), and at max-imum values of the stem cell proliferation rate (p1) andthe reversion rate (r). Increasing the reversion not onlydirectly increases stem cells, but also decreases the level ofHall et al. Journal of Biological Engineering (2017) 11:38 Page 14 of 200 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2Time (min) 104010002000300040005000Cell NumberS P D T0 1 2 3 4 5 6 7 8 9 10Time (min) 1040100020003000400050006000Cell NumberS P D T0 0.5 1 1.5 2 2.5Time (min) 1050123456789Cell Number104S P DFig. 8 Numerical simulations of population dynamics. Top: Population dynamics with initial population of 5000 progenitor cells for C = 7.7 (solid)and C = 6 (dash) with O = 21 andW = 5. Experimental parameters (×10−5): α¯ = 0.01, β¯ = 0.16, γ¯ = 0.16, p¯1 = 119, p¯2 = 160, d¯1 = 13.5, d¯2 = 7.745,r¯ = 0.1.Middle: Population dynamics with initial population of 5000 stem cells and C = 10 for O = 21 (solid) and O = 5 (dash). Same experimentalparameters used as for top simulation. Bottom: Population dynamics after switching to parameters for maximizing D∗ from the initial point of S∗ ,with independent d1, d2, and r. See Fig. 7 for parameters used in maximizing S∗ , and simulation at top for D∗ maximizing parametersnegative feedback to p1 by decreasing the P andD popula-tions, and thus also indirectly contributes to increasing S∗.Similarly, minimizing d1 and d2, the differentiation ratesof stem and progenitor cells, decreases the rate of transferfrom S to P and D, which decreases the level of inhibitoryfeedback to p1 from P and D, indirectly increasing S∗.The progenitor cell population plays a dual role inrelation to the stem cells. It serves as a source of newstem cells via reversion, but inhibits proliferation of stemcells. Thus, it is not clear a priori whether it is better tochoose parameters that lead to a higher or lower equilib-rium progenitor cell population. Above, we chose r andd1 in a way that detracts from the progenitor cell pop-ulation but supports the stem cell population directly.However, the model suggests that p2, the proliferationrate of progenitor cells, and β , the death rate of progen-itor cells, should be chosen to increase the progenitorcell population. In our model, the reversion overcomesthe feedback. However, if the feedback configurationwere modified so that the feedback was stronger, thenHall et al. Journal of Biological Engineering (2017) 11:38 Page 15 of 20maximizing p2 andminimizing β might not be optimal formaximizing S∗. The feedback incorporated into thismodel is only a possible mechanism, not confirmed inbiological reality, so it may be of interest to determineexperimentally whether increasing the proliferation rateand decreasing the death rate of P increases S∗ in reality,which would be consistent with the feedback configu-ration used in this model, or whether it decreases S∗,implying that another feedback mechanism is involved. Ifthe objective is only to maximize the stem cell population,then it would be optimal to maximize γ , the death rateof differentiated cells, which only feed back negatively onproliferation of stem and progenitor cells.In experiments, if we only have control of oxygen level,waste level and contact rate (O,W and C) to alter param-eters, then we cannot achieve all the above goals at once.In particular, γ , like α and β is minimized for oxy-gen levels near 4 or 5%, so if we choose O, W and Cto minimize α and β , then we get a low value of γas well (see Table 2 and Fig. 4). Also, d1 and d2 aremaximized at oxygen levels near 5% but these are notmuch lower at 21% oxygen, whereas other parametersvary more with oxygen level, especially the proliferationrates. So, in the numerical optimization with respect toO in particular, the disadvantage of having a lower deathrate for differentiated cells and slightly higher differen-tiation rates is outweighed by the advantage of havingother parameters at optimal values when oxygen is keptnear 5%.Maximizing P∗ and D∗If our objective is to maximize the progenitor cell popu-lation or the differentiated cell population, it still makessense to keep a strong pool of stem cells, which are asource of progenitors and thus, indirectly, a source of dif-ferentiated cells, and the model confirms that it is best tomaximize the proliferation rates p1 and p2 and minimizethe death rates α and β . If we were only concerned aboutmaximizing P∗, then γ should be maximized again, if wecould do so in practice, and the reversion rate r should beminimized.An initial intuition suggests that maximizing the rateof differentiation from stem cells to progenitor cells (d1)would maximize P∗ because d1 directly transfers cellsfrom S to P. However, because stem cells have a longer lifespan than progenitor and differentiated cells, as reflectedin the values for α, β , and γ , keeping some cells in the stemcell compartment allows for a pool of cells that replenishesthe progenitor cell populations. Thus, the optimal value ofd1 is not necessarily its maximum. Here, the optimal d1for maximizing P∗ and D∗ is an intermediate value withinthe range for d1 (see Table 1).As for d1, the differentiation rate of progenitor cells(d2) is similarly affected by the range in life spans andwe once again get an optimal value of d2 is not necessar-ily it’s maximum. The ability of stem cells to proliferateis also a factor in this result. The lack of proliferationin the differentiated cell population means that it is aterminal state, so driving the system to this compart-ment would drive the populations to extinction eventually.Maintaining S∗ and P∗ populations allows for more Dto be formed without forcing the system to the terminalstate. Thus, the optimal value of d2 is not necessarilyits minimum. Here, the optimal d2 for maximizingD∗ is an intermediate value within the range for d2(see Table 1).Maximizing T∗The total cell population was optimized along with thestem cell population because optimizing the stem cellpopulation decreases the negative feedback on prolifera-tion and, when α is smaller than both β and γ , minimizesthe death rate in the system.Controlling differentiation experimentallyIf differentiation can be controlled chemically in exper-iment, it is no longer dependent on oxygen, waste andcontact rate like the other parameters. In this case, it ispossible to obtain considerably higher steady state popu-lations, whichever subpopulation we wish to optimize (seeindependent results in Table 6).Variable oxygen and wasteThe optimal oxygen concentration in the air above themedium is close to or even lower than physiological levels,depending on which population one intends to maximize,while for waste it is around 2.5 − 5%.When the model parameters are set to match the condi-tions of the experiments, the population sizes and effectsof the change in porosity from 23 to 40% are not quan-titatively captured, as shown by comparing Table 4 andFig. 8. The simulated population size is smaller than thesizes seen experimentally. Also, the large difference inpopulation size noted in experiments is not replicatedby these simulations. However, the simulations do cap-ture the qualitative effect of increasing population size bydecreasing porosity from 23 to 40%. It is likely that thelack of quantitative correlation between the experimentaldata and themodel simulations is a result of the parameterestimations used in the model. This discrepancy in quan-titative results may not be surprising because the dataused to determine both the experimental parameter ratesand the functional effects were very limited. Quantita-tive observations could be captured more effectively withfurther data collection.ConclusionsThe simulated dynamics of our model successfully repro-duce many experimental observations made about howHall et al. Journal of Biological Engineering (2017) 11:38 Page 16 of 20pluripotent stem cells behave when seeded on biomaterialscaffolds. This model allows investigation of alterations inexperimental parameters that would be difficult and costlyto explore in the lab. Analysis of the model confirms theexistence of a unique and stable positive steady state forplausible parameter ranges. A physiological oxygen level isshown to be optimal whichever subpopulation one wishesto maximize. Discrepancies between quantitative aspectsof the model results and experiment could likely be reme-died by improved data on rates of cell behavior and theirfunctional dependence on oxygen, waste and contact.AppendixFlux termsThe explicit flux terms are given byOflux =JO2 · Area of aggregate=[(3.0 × 10−5) 760100 (Oair − O)(5.1737 × 10−8)(60)d]πR2 ,(6)Wflux =JCO2 · Area of aggregate=[(2.5 × 10−5) 760100 (Wair − W )(5.1737 × 10−8)(60)d]πR2(7)where d = H − R is the distance from the top of themedium to the center of the aggregate,H = Vmedium +43πR39.6is the height of the medium,R =( S + P + D8 · 4500 · 1000)1/3is the radius of the aggregate, and where the base of theculture plate has an area of 9.6cm2. Note that 3.0 × 10−5and 2.5 × 10−5 are the diffusion coefficients for O andW,respectively. The remaining constants in the numeratorsof the flux terms are unit conversion terms.Cell body cluster area to cell number calculationThe number of cells was calculated from cell body clus-ter area assuming an initial population of 4500 cells.The aggrewell platform used produces aggregates withbetween 4000 and 5000 cells. With this assumption inplace, the calculation of cell numbers from cell bodycluster area was determined as follows.For any time t, we have that the approximate number ofcells in the aggregate, N(t), is given byN(t) = CDo · V (t) =(6(4500)πcellsmm3)(43π√A(t)π)where A(t) is the experimentally determined cell bodycluster area, V (t) is the volume of the aggregate (derivedfrom A(t)), and CDo is the cell density of the aggregate,which is assumed to be constant. We first determine theinitial volume of cells using an averaged cell body clusterarea of 0.8mm2:0.8mm2 = A(0) = πr(0)2 =⇒ r(0) = 0.5mm=⇒ V (0) = 43π(0.5mm)3 ,where r(0) is the radius of the aggregate at t = 0. Thus,the cell density is given byCDo =4500 cellsV (0) =6(4500) cellsπmm3 .Finally, to convert cell body cluster area to number ofcells, we have N(t) = CDo · V (t) .Themodel without stem cellsThis version of the model corresponds to the larger modelin the case where the initial stem cell population is zeroand the reversion rate, r, is zero. This reduces it to a two-variablemodel, for whichmore analytic tools are available.The simplified PDmodel isdPdt = − βP +2p2P1 + P + D − d2PdDdt = − γD + d2P .(8)It is more analytically tractable than the SPD system, butour analysis strategy is the same.First, it can be shown thatP ∈[0, 2p2β + d2]and D ∈[0, 2p2d2γ (β + d2)](9)is invariant for the PD system and that solutions start-ing outside this region in the non-negative quadrant fallinto it. Also, there are no periodic solutions by Dulac’sCriterion [36].Existence and stability of fixed pointsThe nullclines of P and D, the intersection of which deter-mine the fixed points, are given by setting the right handsides of Eq. (8) to zero. Solving this system gives two fixedpoints: (0, 0), and a unique positive fixed point (P∗,D∗),which exists if and only if 2p2 > β + d2, whereP∗ = γ (2p2 − β − d2)(β + d2)(γ + d2) , (10)D∗ =d2(2p2 − β − d2)(β + d2)(γ + d2) , (11)T∗ =P∗ + D∗ = 2p2β + d2 − 1 . (12)The existence condition for the positive fixed point isfeasible within the parameter space defined in Table 2.Hall et al. Journal of Biological Engineering (2017) 11:38 Page 17 of 20Also, in our parameter space min p2 = 0, which meansthat the situation where the positive fixed point does notexist is also possible.The Jacobian matrix of the PD system isJ(P,D) =[−β − d2 + 2p2(1+D)(1+P+D)2 − 2p2P(1+P+D)2d2 −γ]and it is easy to check that the zero equilibrium is a sta-ble node for 2p2 < β + d2 and is a saddle point, andtherefore unstable, for 2p2 > β + d2. It is also easy toshow that Trace(J) < 0 and Det(J) > 0 for the pos-itive equilibrium (P∗,D∗) whenever it exists, i.e., when2p2 > β + d2. In fact, the positive fixed point is globallystable with respect to the invariant region (9) apart fromthe D-axis. This follows from the lack of periodic solu-tions and the Poincaré-Bendixson Theorem applied to theinterior of the invariant region. Thus, when the positivefixed point exists, all solutions converge to it. Otherwise,all solutions go to zero.It should be noted that although our model precludessolutions reaching the D-axis from elsewhere when thepositive fixed point exists (i.e., precludes an initially posi-tive P from reaching 0), the model is not really appropriatefor very small populations, so if P = 1 or a very smallnumber, then in reality we might expect stochastic effectsto allow P to jump to 0. This is not a real problem, how-ever, for the ranges of cell numbers we deal with in ourexperimental setup.Optimizing the positive PD fixed pointAs for the SPD model, we consider each subpopulation,or the total population, at the fixed point as a functionof each parameter, and explore the effects of individualparameters on the steady state (sub)populations.From Eq. (12), the total cell population at equilibrium,T∗, is maximized when p2 is maximized, and when β andd2 are minimized. The parameter γ has no effect on T∗.The progenitor cell population at equilibrium, P∗, givenby Eq. (10) is increasing in p2 and γ , and decreasing in βand d2, as can be determined by considering the sign of thederivative with respect to each parameter. Thus, the pro-genitor cell population at equilibrium, P∗, is maximizedwhen p2 and γ are maximized, and when β and d2 areminimized.Equation (11) shows that the differentiated cell popula-tion at equilibrium, D∗, is increasing in p2 and decreasingin β and γ , but is neither a strictly increasing or decreas-ing function of d2. Solving dD∗dd2 = 0 for positive values ofd2 givesd2 = −βγ +√β2γ 2 + βγ (2p2 + γ )(2p2 − β)2p2 + γat which D∗ has a maximum. We denote this value as d∗2,which, for most values of p2, γ , and β , lies in the allowablerange for d2.The optimal values of the parameters for the PD modelare given in Table 7. The optimized fixed points, that is,the values calculated by including the optimal parameterswithin our given parameter ranges in Eqs. (10) and (11),are given in Table 8. The populations under optimizationof T∗ are given both for the minimum and maximum val-ues of γ , since T∗ does not depend on γ , though P∗ andD∗ do.The above fixed points are not experimentally realizablesince the parameters are not independent, so we appliedthe same optimization process as for the SPD system tak-ing into account the dependence of the parameters on O,W, and C (Table 9).Including variable oxygen and waste in the modelwithout stem cells The PDOW system is given bydPdT = − βP +2p2P1 + P + D − d2PdDdt = − γD + d2PdOdt = − u2P − u3D + Oflux(P,D,O)dWdt =w2P + w3D + Wflux(P,D,W ) .(13)In the PD system, we determined the optimal popula-tions and the associated parameters, as well as optimaloxygen and waste concentrations. Using the PD fixedpoints and including the optimal O and W values as theO andW components of the PDOW fixed point, the fixedpoint equations for O and W can be solved for Oair andWair .The fixed points for P and D in the PDOW systemare still those given in Eqs. (10) and (11), and the cor-responding value of T∗ is still the γ -independent valuegiven in Eq. (12). Numerical values of all the variablesat fixed points corresponding to the various optimizationschemes, are given in Table 9. It can also be shown (fromeigenvalues of the Jacobian matrix) that these fixed pointsare all locally asymptotically stable. The Oair and WairTable 7 Summary of optimal parameter values for maximizingcell populationsParameter Max P∗ Max D∗ Max T∗β Min Min Minγ Max Min No effectp2 Max Max Maxd2 Min d∗2 MinHall et al. Journal of Biological Engineering (2017) 11:38 Page 18 of 20Table 8 Optimized fixed points with O-independent parametersOptimized value P∗ D∗ T∗ D∗/T∗P∗ 3325.2 11.7 3336.9 0.0035D∗ 857.5 993.3 1850.8 0.5367T∗ , min γ 3194.2 142.7 3336.9 0.0428T∗ , max γ 3325.2 11.7 3336.9 0.0035values required to achieve the optimal O and W at theequilibrium values of P and D are also given in Table 9.Discussion of the model without stem cellsFixed point existence and stability The model withoutstem cells always has a fixed point at zero, as expected ina population model. When 2p2 < β + d2, the zero fixedpoint is stable and there are no other equilibria in the fea-sible region, so the system will eventually decay to 0 fromany initial point. Biologically, 2p2 ≤ β+d2 corresponds tothe situation where the production of progenitor cells byproliferation is overcome by the loss of progenitor cells bydeath and differentiation. Differentiated cells are a termi-nal state and inhibit the proliferation of progenitor cells,so the system is driven to extinction.At 2p2 = β + d2, there is a bifurcation above whichthe zero fixed point becomes unstable and a stable pos-itive fixed point, (P∗,D∗), is formed. Thus, for 2p2 >β + d2, both progenitor and differentiated cells persist.In this case, the proliferation rate overcomes the loss ofprogenitor cells by death and differentiation, allowing forthe maintenance of a pool of progenitor cells. This poolof progenitor cells also maintains a population of dif-ferentiated cells as the progenitor cells transition to thedifferentiated state.Fixed point dependence on parameters The progenitorcell population at equilibrium, P∗, is maximized when βand d2, the death rate and differentiation rate of progeni-tors, are minimized, and when p2 and γ , the proliferationrate of progenitors and the death rate of differentiatedcells, are maximized. The total cell population, T∗, is max-imized under the same conditions except that in this casethe death rate of differentiated cells has no effect. Thedifferentiated cell population at equilibrium, D∗, is maxi-mized when β and γ , the two death rates, are minimized,when p2, proliferation of progenitors, is maximized, andwhen the differentiation rate is at an optimal intermediatevalue, d2 = d∗2. The proportion of D∗ is maximized whenγ is minimized and d2 is maximized, while p2 and β donot affect the proportion.Progenitor cells are the source of new cells in the 2-dimensional system, so the optimal value of p2 and βfor optimization of all populations are the maximum andminimum values, respectively. The optimal values of p2and β have less complicated implications than the otherparameters as progenitor cells only feed back on theirown proliferation. To maximize any cell population, max-imizing the production and minimizing the removal ofthe proliferative progenitor cell population must be opti-mal as it supplies both progenitor and differentiated cells.Decreasing the differentiation rate of progenitor cells willincrease both progenitor and total cell populations. Thiseffect is a result of two factors. First, the differentiatedcell state is terminal. Therefore, prolonging the length oftime the cells stay in the proliferative progenitor state (bydecreasing d2) increases the rate of production of newcells. Secondly, the differentiated cell population nega-tively feeds back onto the proliferation of the progeni-tor cells. Thus, as d2 increases and progenitor cells aredriven to the differentiated cell population, the inhibi-tion of progenitor cell proliferation increases and will leadto a smaller population of progenitor cells and total cellpopulation.A different d2 value maximizes the differentiated cellpopulation: d2 = 0.0000094522. This numerical valuedepends on the values of the other parameters. Thus, if theexperimental ranges for parameters were altered so thatTable 9 Optimized fixed points for the PD system with dependent parametersOptimized value P∗ D∗ T∗ D∗/T∗ O∗ Oair WairP∗ , coupled 141.6 123.0 264.6 0.4648 0.8976 1.2813 4.5312P∗ , d2 uncoupled 2856.1 53.2 2909.3 0.0183 4.2896 4.9931 4.1407D∗ , coupled 53.3 330.7 384.1 0.8611 4.0059 4.4573 4.4485D∗ , d2 uncoupled 908.2 713.2 1621.4 0.4399 4.3585 4.9996 4.2168T∗ , coupled, γmin 60.4 345.0 405.3 0.8511 3.7393 4.1973 4.4404T∗ , coupled, γmax 89.7 315.6 405.3 0.7786 3.7393 4.1932 4.4454T∗ , d2 uncoupled, γmin 2823.8 85.5 2909.4 0.0294 4.2891 4.9973 4.1394T∗ , d2 uncoupled, γmax 2856.1 53.2 2909.4 0.0183 3.7393 4.9926 4.1407Note thatW∗ = 5 in all cases. “Coupled” indicates dependence of the parameters on O,W and C. “Uncoupled” indicates that the independent optimal parameter value wasused. For the optimal T∗ populations, the populations are given for the minimum and maximum values of γ , γmin and γmax, because T∗ does not depend on γ , though P∗and D∗ doHall et al. Journal of Biological Engineering (2017) 11:38 Page 19 of 20the optimal values of the parameters were also changed,this optimal d2 value would be different. For our set ofexperimental parameters, the optimal d2 values is lowwithin the range for the compound d2 parameter. Thisresult may not be immediately obvious because an an ini-tial intuition might be to maximize the rate of creation ofdifferentiated cells in order to maximize D∗. However, thedifferentiated cell state is a terminal state, so the optimald2 occurs where there is a balance between creation ofdifferentiated cells and the maintenance of the progenitorcell pool that leads to differentiated cells. The differen-tiated cell population is maximized by keeping a sizableprogenitor pool to produce more differentiated cells. Thisprogenitor cell population ismaintained by not driving thepopulation too quickly to the differentiated state as wellas keeping the progenitor cell population large enough toovercome the proliferative inhibition from the differenti-ated cell population. Thus, the size of D∗ is maximized atthis intermediate value of d2.The parameter γ is also somewhat counterintuitive.Increasing the death rate of differentiated cells increasesthe progenitor cell population because fewer differen-tiated cells mean less inhibition of proliferation of theprogenitor cells. Minimizing γ also maximizes D∗. Thisis not as obvious as it seems because of the feedback.Although a lower death rate increases the differentiatedcell population, it also increases the negative feedbackon the proliferation of progenitor cells, which are thesource of differentiated cells. Thus, minimizing γ is nota priori the correct method for maximizing D∗. How-ever, with the current feedback mechanism, the decreasedloss of differentiated cells outweighs the decrease in inputfrom the progenitor cell population. Therefore, at leastwith the current feedback mechanism, minimizing γmaximizes D∗.The results of optimization for p2 and β are obviousbiologically and are not complicated by other factors: toincrease a population, increase proliferation and decreasedeath. The model is consistent with this observation. Themore counterintuitive results for d2 and γ provide moreinsight into the dynamics of the system.Fixed point optimization The optimization of theparameters was first carried out assuming each parame-ter could be altered independently, in order to understandhow the populations were affected by each parameter.However, the parameters are actually determined by oxy-gen, waste and scaffold porosity, and as such cannot beconsidered independent. Thus, to provide useful informa-tion to experimentalists, the populations were optimizedwith respect to O,W and C by determining which val-ues resulted in the closest parameter set to the optimalset. The resulting O values indicated that a physiolog-ical O2 level, or lower, is the best for maximizing thepopulations because most parameters were optimized atthese O values, while d2 was the only clear exception. Asd2 was an outlier, the above optimization was also carriedout when d2 was optimized as an independent parame-ter (corresponding to altering the differentiation rate ofprogenitors by appropriate chemical factors in experi-ments). The resulting optimal oxygen level was higherwhen d2 was uncoupled, although still within the physio-logical range. It is also of note that the population levelswith d2 decoupled are significantly higher than the fullycoupled counterparts. The optimal W was determined tobe lower than the typical CO2 levels used experimentally.The contact rate C determined to be optimal for all popu-lations was 3.75, which agrees in a general sense with thefit to the experimental data.The population sizes of the fixed points are muchsmaller than the populations in the experimental data.This may be a consequence of the limited data used forboth rate determination and fitting, and consequent inac-curacy in those values. Another possible reason for thisdiscrepancy is that stem cells were actually present in theexperiments, despite the intention of seeding only progen-itors. The presence of stem cells in the seeded populationtends to increase the total cell numbers over time.AbbreviationshiPSC: Human induced pluripotent stem cell; iPSC: Induced pluripotent stemcell; NIM: Neural induction medium; ODE: Ordinary differential equation; PCL:Poly-(-caprolactone)AcknowledgmentsThe authors would like to acknowledge Junghyuk Ko and Martin Byung-GukJun for their contributions to the laboratory work.FundingThe study was supported by Natural Sciences and Engineering ResearchCouncil (NSERC) Discovery Grants.Availability of data andmaterialsPlease contact author for data requests.Author’s contributionsMH and RE responsible for mathematical model development and analysis.MH responsible for numerical simulations. NKM responsible for all laboratorywork. MH, RE, and SW participated in drafting the manuscript. All authors readand approved manuscript.Ethics approval and consent to participateAll appropriate ethical approvals were obtained for conducting this research.Consent for publicationAll authors have approved the manuscript and agreed with its submission.Competing interestsThe authors declare that they have no competing interests.Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims inpublished maps and institutional affiliations.Author details1Department of Mathematics and Statistics, University of Victoria, Victoria,Canada. 2Department of Mechanical Engineering, University of Victoria,Victoria, Canada. 3Division of Medical Sciences, University of Victoria, Victoria,Canada. 4Department of Biochemistry, University of British Columbia,Hall et al. Journal of Biological Engineering (2017) 11:38 Page 20 of 20Vancouver, Canada. 5International Collaboration on Repair Discoveries,University of British Columbia, Vancouver, Canada. 6Centre for BiomedicalResearch, University of Victoria, Victoria, Canada.Received: 22 May 2017 Accepted: 20 September 2017References1. Langer R, Vacanti J. Tissue engineering. Science. 1993;260(5110):920–6.2. Willerth SM, Sakiyama-Elbert SE. Approaches to neural tissue engineeringusing scaffolds for drug delivery. Adv Drug Deliv Rev. 2007;59(4):325–38.3. Narsinh KH, Plews J, Wu JC. Comparison of human induced pluripotentand embryonic stem cells: Fraternal or identical twins? Mol Ther.2011;19(4):635–8.4. Takahashi K, Yamanaka S. Induction of pluripotent stem cells from mouseembryonic and adult fibroblast cultures by defined factors. Cell.2006;126(4):663–76.5. Keung AJ, de JuanPardo EM, Schaffer DV, Kumar S. Rho gtpases mediatethe mechanosensitive lineage commitment of neural stem cells. StemCells. 2011;29(11):1886–97.6. Keung AJ, Dong M, Schaffer DV, Kumar S. Pan-neuronal maturation butnot neuronal subtype differentiation of adult neural stem cells ismechanosensitive. Sci Rep. 2013;3:1817.7. Roccio M, Schmitter D, Knobloch M, Okawa Y, Sage D, Lutolf MP.Predicting stem cell fate changes by differential cell cycle progressionpatterns. Development (Cambridge, England). 2013;140(2):459–70.8. Studer L, Csete M, Lee SH, Kabbani N, Walikonis J, Wold B, McKay R.Enhanced proliferation, survival, and dopaminergic differentiation of cnsprecursors in lowered oxygen. J Neurosci. 2000;20(19):7377–83.9. Lane SW, Williams DA, Watt FM. Modulating the stem cell niche fortissue regeneration. Nat Biotechnol. 2014;32(8):795–803.10. Lo WC, Chou CS, Gokoffski KK, Wan FY-M, Lander AD, Calof AL, Nie Q.Feedback regulation in multistage cell lineages. Math Biosci Eng.2009;6(1):59–82.11. Herberg M, Roeder I. Computational modelling of embryonic stem-cellfate control. Development. 2015;142(13):2250–260.12. Jin K, Mao XO, Batteur S, Sun Y, Greenberg DA. Induction of neuronalmarkers in bone marrow cells: differential effects of growth factors andpatterns of intracellular expression. Exp Neurol. 2003;184(1):78–89.13. Willerth SM. Neural tissue engineering using embryonic and inducedpluripotent stem cells. Stem Cell Research & Therapy. 2011;2:17.14. Mohtaram NK, Ko J, King C, Sun L, Muller N, Jun MB-G, Willerth SM.Electrospun biomaterial scaffolds with varied topographies for neuronaldifferentiation of human-induced pluripotent stem cells. J Biomed MaterRes A. 2015;103(8):2591–601.15. Takahashi K, Yamanaka S. Induction of pluripotent stem cells from mouseembryonic and adult fibroblast cultures by defined factors. Cell.2006;126(4):663–76.16. Pagliara S, Franze K, McClain CR, Wylde GW, Fisher CL, Franklin RJM,Kabla AJ, Keyser UF, Chalut KJ. Auxetic nuclei in embryonic stem cellsexiting pluripotency. Nat Mater. 2014;13(6):638–44.17. Robinson M, Yau S-y, Sun L, Gabers N, Bibault E, Christie BR, Willerth SM.Optimizing differentiation protocols for producing dopaminergic neuronsfrom human induced pluripotent stem cells for tissue engineering applications.Biomarker Insights. 2015;10(S1):61–70.18. Ko J, Mohtaram NK, Ahmed F, Montgomery A, Carlson M, Lee PC,Willerth SM, Jun MB. Fabrication of poly (-caprolactone) microfiberscaffolds with varying topography and mechanical properties for stemcell-based tissue engineering applications. J Biomater Sci Polym Ed.2014;25(1):1–17.19. Engler AJ, Sen S, Sweeney HL, Discher DE. Matrix elasticity directs stemcell lineage specification. Cell. 2006;126(4):677–89.20. Mohyeldin A, Garón-Muvdi T, Quiõnes-Hinojosa A. Oxygen in stem cellbiology: A critical component of the stem cell niche. Cell Stem Cell.2010;7(2):150–61.21. Vieira HLA, Alves PM, Vercelli A. Modulation of neuronal stem celldifferentiation by hypoxia and reactive oxygen species. Prog Neurobiol.2011;93(3):444–55.22. Zhu J, Aja S, Kim E, Park MJ, Ramamurthy S, Jia J, Hu X, Geng P,Ronnett GV. Physiological oxygen level is critical for modeling neuronalmetabolism in vitro. J Neurosci Res. 2012;90(2):422–34.23. Braam SR, Zeinstra L, Litjens S, Ward-van Oostwaard D, van den Brink S,van Laake L, Lebrin F, Kats P, Hochstenbach R, Passier R, et al.Recombinant vitronectin is a functionally defined substrate that supportshuman embryonic stem cell self-renewal via αvβ5 integrin. Stem Cells.2008;26(9):2257–265.24. Montgomery A, Agbay A, Edgar J, Gabers N, Gomez J, Mohtaram N,King C, Mitchell A, Rajwani A, Rattray D, Robinson M, Shapka A, Sun L,Wong A, Willerth S. Combining protein-based biomaterials with stemcells for spinal cord injury repair. OA Stem Cells. 2014;2(1):1.25. Montgomery A, Wong A, Gabers N, Willerth SM. Engineeringpersonalized neural tissue by combining induced pluripotent stem cellswith fibrin scaffolds. Biomater Sci. 2015;3(2):401–13.26. Clarke L, van der Kooy D. Low oxygen enhances primitive and definitiveneural stem cell colony formation by inhibiting distinct cell deathpathways. Stem Cells (Dayton, Ohio). 2009;27(8):1879–86.27. Ezashi T, Das P, Roberts RM. Low o2 tensions and the prevention ofdifferentiation of hes cells. Proc Natl Acad Sci U S A. 2005;102(13):4783–788.28. Gray DR, Chen S, Howarth W, Inlow D, Maiorella BL. Co2 in large-scaleand high-density cho cell perfusion culture. Cytotechnology.1996;22(1-3):65–78.29. Birket MJ, Orr AL, Gerencser AA, Madden DT, Vitelli C, Swistowski A,Brand MD, Zeng X. A reduction in atp demand and mitochondrial activitywith neural differentiation of human embryonic stem cells. J Cell Sci.2011;124(Pt 3):348–58.30. Cho CH, Park J, Nagrath D, Tilles AW, Berthiaume F, Toner M, Yarmush ML.Oxygen uptake rates and liver-specific functions of hepatocyte and 3t3fibroblast co-cultures. Biotech Bioeng. 2007;97(1):188–99.31. Salomoni P, Calegari F. Cell cycle control of mammalian neural stem cells:putting a speed limit on g1. Trends Cell Biol. 2010;20(5):233–43.32. Hardwick LJ, Philpott A. Nervous decision-making: to divide ordifferentiate. Trends Genet. 2014;30(6):254–61.33. Altman BJ, Rathmell JC. Metabolic stress in autophagy and cell deathpathways. Cold Spring Harb Perspect Biol. 2012;4(9):008763.34. Vander Heiden MG, Cantley LC, Thompson CB. Understanding thewarburg effect: the metabolic requirements of cell proliferation. Science.2009;324(5930):1029–33.35. Scientific TF. pH and CO2 Levels. https://www.thermofisher.com/ca/en/home/references/gibco-cell-culture-basics/cell-culture-environment/ph-co2-levels.html. Accessed 29 June 2016.36. Strogatz SH. Nonlinear Dynamics and Chaos: with Applications to Physics,Biology, Chemistry, and Engineering. Boulder, CO: Westview press; 2014.37. Guarino RD, Dike LE, Haq TA, Rowley JA, Pitner JB, Timmins MR. Methodfor determining oxygen consumption rates of static cultures frommicroplate measurements of pericellular dissolved oxygen concentration.Biotechnol Bioeng. 2004;86:775–87.• We accept pre-submission inquiries • Our selector tool helps you to find the most relevant journal• We provide round the clock customer support • Convenient online submission• Thorough peer review• Inclusion in PubMed and all major indexing services • Maximum visibility for your researchSubmit your manuscript atwww.biomedcentral.com/submitSubmit your next manuscript to BioMed Central and we will help you at every step:
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Faculty Research and Publications /
- Modeling the behavior of human induced pluripotent...
Open Collections
UBC Faculty Research and Publications
Modeling the behavior of human induced pluripotent stem cells seeded on melt electrospun scaffolds Hall, Meghan E; Mohtaram, Nima K; Willerth, Stephanie M; Edwards, Roderick Oct 23, 2017
pdf
Page Metadata
Item Metadata
Title | Modeling the behavior of human induced pluripotent stem cells seeded on melt electrospun scaffolds |
Creator |
Hall, Meghan E Mohtaram, Nima K Willerth, Stephanie M Edwards, Roderick |
Contributor | International Collaboration on Repair Discoveries |
Publisher | BioMed Central |
Date Issued | 2017-10-23 |
Description | Background: Human induced pluripotent stem cells (hiPSCs) can form any tissue found in the body, making them attractive for regenerative medicine applications. Seeding hiPSC aggregates into biomaterial scaffolds can control their differentiation into specific tissue types. Here we develop and analyze a mathematical model of hiPSC aggregate behavior when seeded on melt electrospun scaffolds with defined topography. Results: We used ordinary differential equations to model the different cellular populations (stem, progenitor, differentiated) present in our scaffolds based on experimental results and published literature. Our model successfully captures qualitative features of the cellular dynamics observed experimentally. We determined the optimal parameter sets to maximize specific cellular populations experimentally, showing that a physiologic oxygen level (∼ 5%) increases the number of neural progenitors and differentiated neurons compared to atmospheric oxygen levels (∼ 21%) and a scaffold porosity of ∼ 63% maximizes aggregate size. Conclusions: Our mathematical model determined the key factors controlling hiPSC behavior on melt electrospun scaffolds, enabling optimization of experimental parameters. |
Subject |
Stem cell Ordinary differential equation Differentiation Proliferation Tissue engineering |
Genre |
Article |
Type |
Text |
Language | eng |
Date Available | 2018-05-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution 4.0 International (CC BY 4.0) |
DOI | 10.14288/1.0366880 |
URI | http://hdl.handle.net/2429/65933 |
Affiliation |
Medicine, Faculty of Other UBC Non UBC Biochemistry and Molecular Biology, Department of |
Citation | Journal of Biological Engineering. 2017 Oct 23;11(1):38 |
Peer Review Status | Reviewed |
Scholarly Level | Faculty |
Copyright Holder | The Author(s) |
Rights URI | http://creativecommons.org/licenses/by/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 52383-13036_2017_Article_80.pdf [ 2.04MB ]
- Metadata
- JSON: 52383-1.0366880.json
- JSON-LD: 52383-1.0366880-ld.json
- RDF/XML (Pretty): 52383-1.0366880-rdf.xml
- RDF/JSON: 52383-1.0366880-rdf.json
- Turtle: 52383-1.0366880-turtle.txt
- N-Triples: 52383-1.0366880-rdf-ntriples.txt
- Original Record: 52383-1.0366880-source.json
- Full Text
- 52383-1.0366880-fulltext.txt
- Citation
- 52383-1.0366880.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52383.1-0366880/manifest