{"paper_id":"221d63aa-cd00-4b64-b3ee-bef86334ae58","body_text":"Mechanistic interpretation of biological tissue\ngrowth experiments with a computational model\nShahak Kuba1, Matthew J. Simpson1,2, Pascal R. Buenzli1,2∗\n1School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia.\n2ARC Centre of Excellence for the Mathematical Analysis of Cellular Systems, QUT, Brisbane,\nAustralia.\nMarch 13, 2026\nAbstract –The growth rates of biological tissues are influenced by the existing substrate\ngeometry, mechanobiological processes and the interplay between them. Disentangling the\ncontributions of geometry and mechanobiology experimentally is challenging, as mechanical\nproperties are difficult to measure and tissue samples provide only static snapshot in time.\nHowever, the composition of a tissue preserves cues of the dynamic processes that shaped its\narchitecture. In this work, we present a computational model of tissue growth that captures\naspects of the interplay between geometry, mechanics, and stochastic biological processes,\nwhich we use to generate synthetic tissue compositions directly comparable with experimental\nsamples. This framework enables quantitative analysis of tissue morphology, inference of\nunderlying growth mechanisms, estimation of dynamic rates from single-time-point data, and\ninvestigation of how stochasticity contributes to emergent growth patterns. We demonstrate\nthe applicability of the model to simulate the growth of different tissue types by applying this\nframework to two distinct tissue growth scenarios: (i) tissue grown within 3D-printed porous\nscaffolds, and (ii) bone formation in cortical pores.\nKey words:Tissue Engineering, cell culture, bone, single-cell, stochastic\n1 Introduction\nThegrowthofbiologicaltissuesiscentraltothedevelopment,repair,andregenerationoforgans\nand tissues. Tissue growth plays an important role for example in bone remodeling, tumor\ngrowth,andwoundhealing(Martinetal., 1998;Poujadeetal., 2007;Rollietal., 2012;Wozniak\nand El Haj,2007; Lowengrub et al.,2009), as well as in tissue engineering and regenerative\nmedicine,wherebiologicaltissuegrowthisinvolvedinengineeringtherepairorreplacementof\ndamaged or diseased tissues (Hollister,2005;O’Brien,2011;Bidan et al.,2013). During tissue\n∗Corresponding author. Email address:pascal.buenzli@qut.edu.au\n1\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigure 1– Images from case studies of biological tissue growth. (a)In vitrotissue growth of\nMC3T3 cells within a400 μm × 400 μm 3D-printed square porous scaffold cultured over 14\ndays. The extracellular matrix (ECM) is depicted in green and differentiated cell nuclei are\nindicated in blue. (b)In vivocross section of human bone tissue with multiples cortical pores\nhighlightedinredsurroundingtheHaversiancanalsystemwithembeddedcells(osteocytes)are\nindicated in black (reproduced with permission fromSkedros et al.(2005)).\ngrowth, new biological material is usually created within constrained spaces with experimental\nobservations indicating that the geometry of both the constraining space and the existing tissue\nsubstrate influences the rate of growth (Rumpler et al.,2008; Bidan et al.,2012; Schamberger\net al.,2023). The accumulation of new material within these geometric constraints induces\ncrowding effects and mechanical stresses on individual cells, which may affect their ability\nto proliferate, differentiate, and survive through mechanobiological processes (Nelson et al.,\n2005; Discher et al.,2005; Lim et al.,2010; Ladoux and Mége,2017). Additionally, cell-scale\ncues arising from the geometry of the existing substrate are also known to influence cell fate\n(McNamaraetal., 2010;Dobbengaetal., 2016;Callensetal., 2020). Thetissuegrowthdynamics\narising from the interplay between geometry, mechanics, and cell fate in turn influence the rate\noftissuegrowth,yettheprecisequantitativedetailsoftheseinterplayremainpoorlyunderstood.\nExperimental approaches for studying the interactions between geometry and mechanics can\noften be challenging, especially as tissue samples provide single snapshots in time and offer\nlimitedinsightintogrowthdynamics. However,thecompositionofatissuecontainsusefulcues\nthat provide a record of some of the dynamic processes that led to the architecture of a tissue.\nBiologicaltissuesareprimarilycomposedofanextracellularmatrix(ECM)anddifferentiated\ncellsembeddedwithinit(Albertsetal., 2002;KurniawanandBouten, 2018). Sincethematerial\ncompositionformsduringgrowth,thepropertiesofdifferentiatedcellsandECMencodevaluable\ninformation about the developmental history and can provide a fingerprint of the underlying\nprocesses that are regulated in healthy tissues, and dysregulated in diseased tissues. Confocal\nmicroscopy and quantitative backscattered electron imaging of tissue samples, such as those in\nFigure1,provideameanstomeasurematerialproperties,includingspatialinformationrelating\nto differentiated cells. Several studies use these measurements to characterize ECM properties\nsuch as fiber alignment as well as the shape, orientation, and number densities of differentiated\ncellswithindifferenttissuetypes(Roschgeretal., 2008;Vatsaetal., 2008;Ascenzietal., 2008;\nvan Hove et al.,2009; Britz et al.,2012; Mader et al.,2013; Lanaro et al.,2021). Beyond\n2\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nproviding a quantitative description of biological tissue composition, these data sets can be\nleveragedwithinmathematicalmodelingtohelpinferthemechanismsthatgoverntissuegrowth\nsincethesemodelscanbedesignedtoexplicitlyincorporatespecificmechanisms. Forexample,\nmineral density distributions in trabecular bone have been used to study bone mineralization\ndynamics (Ruffoni et al.,2007), osteocyte density in cortical pores to estimate osteoblast burial\nrates(Buenzli, 2015), andmouseincisorenamelmeasurementstoformulatepredictivetheories\nofcellularresponsestostrainvariations(Coxetal., 2024). Inthisregard,mathematicalmodeling\noffersapowerfulframeworkforgeneratingandtestinghypotheses,analyzingexperimentaldata,\nand could be useful for disentangling interplay between geometry, mechanics, and cell fate in\nthe processes underlying tissue growth.\nSeveral mathematical models of tissue growth have been developed to predict and help\nunderstand the evolution of biological tissues, with varying levels of biological detail. Phe-\nnomenologicalmodelscapturethemacroscopicbehaviorofgrowingtissuesbydefininggrowth\nlawsthatlinkthevelocityofthetissueboundaryanditscurvature,forexample,meancurvature\nflow (Rumpler et al.,2008; Guyot et al.,2014; Callens et al.,2023). Continuum mechan-\nics approaches typically employ growth tensors to describe tissue deformation and expansion\n(Ambrosi and Guana,2007; Dunlop et al.,2010; Goriely, 2007). However, these growth laws\nand growth tensors are difficult to relate to biological processes. Cell-based models, which\ninclude both continuum and discrete approaches, complement phenomenological descriptions\nby explicitly accounting for cell populations and modeling their associated behaviors (Bidan\net al.,2012; Joly et al.,2013; Gahffarizadeh et al.,2018). These modeling frameworks have\nshown that some curvature effects arise from the crowding and spreading of cells (Alias and\nBuenzli,2017;Hegarty-Cremeretal., 2021;Kubaetal., 2026). Attheindividualcelllevel,these\nmodels are also able to connect mechanical stresses to tissue organization and to cellular events\n(Baker et al.,2019; Murphy et al.,2019; Buenzli et al.,2025; Brown et al.,2025). While these\nmodeling approaches have provided valuable insights into tissue growth, they tend to focus on\ntheevolutionofthetissuesurfaceratherthangeneratingatissuecompositionthatiscomparable\nto experimental samples.\nInthiswork,wepresentacomputationalmodelofbiologicaltissuegrowththatcanbeusedto\nsimulateandanalyzethedynamicprocessesthatleadtoobservedtissuecompositionsinconfined\ngeometries. The model explicitly represents individual cells and includes mechanical interac-\ntions between cells, the generation of ECM, and incorporates cell proliferation, differentiation,\nand death as stochastic events. Unlike previous interface-based or curvature-controlled growth\nmodels, this framework explicitly records the positions and orientations of differentiated cells\nas tissue forms, enabling direct comparison of synthetic and experimental tissue compositions\nat single-cell resolution. Comparisons between model-generated predictions and experimental\ndataenablesusto: (i)relatemechanicalinteractionsbetweencellstoboththeirorientationsand\nthe shape of the tissue interface; (ii) identify key mechanisms that are important for regulating\ntissue growth; (iii) estimate cell differentiation rates at different times from single snapshot;\n3\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nand,(iv)analyzehowmodeledmechanismsandstochasticityshapetissuecompositionandgive\nrise to random growth patterns. To demonstrate the flexibility of this framework, we apply it\ntoinvestigatedatafromtissueengineeringcellcultureexperimentsin3D-printedscaffolds,and\nbone formation in cortical pores.\n2 Materials & Methods\nThis study is primarily simulation-based, using comparisons between computational model\npredictionsandexperimentalobservationstointerpretthemechanismsunderlyingtissuegrowth.\nWhile experimental data provide quantitative description of tissue composition, they do not\nreveal the precise mechanisms that control the growth behaviour. By applying mechanistic\nsimulation models to interpret these data, we obtain biological insights that are otherwise\ndifficult to access.\n2.1 Mathematical model\nThe mathematical model simulates tissue growth by evolving a discrete representation of the\ntissueinterface,withnewlyformedtissuefillingtheregionbehindit. Theinterfaceisrepresented\nas a one-dimensional chain of𝑁 tissue-forming cells connected through their boundaries. The\npositionofthe 𝑖thcellattime 𝑡isgivenbythetwo-dimensionalpositionvectorsofitsboundaries,\n𝒓𝑖−1(𝑡) and 𝒓𝑖 (𝑡) (Figure2a).\nTissuegrowtharisesfromthreeprocesses: (i)mechanicalinteractionsbetweenneighbouring\ncells, (ii) ECM deposition at the interface, and (iii) stochastic cellular events. Following\nprevious discrete cell–cell mechanical interaction modeling approaches, intercellular forces are\nrepresented by representing cell bodies as mechanical springs (Murray et al.,2009; Murphy\net al.,2019; Buenzli et al.,2025). Cells are assumed to generate an area of tissue material\nat a constant rate𝑘f (μm2/day) that displaces them perpendicularly to the interface (Buenzli,\n2015; Alias and Buenzli,2017; Kuba et al.,2026). Together, these two mechanisms yield two\nindependentvelocitycontributionsthatevolvethepositionofagivencellboundary 𝒓𝑖,according\nto\nd𝒓𝑖\nd𝑡 =𝒗 m\n𝑖 +𝒗 f\n𝑖, 𝑖=1,2, ..., 𝑁,(1)\nwhere𝒗 m\n𝑖 represents tangential motion arising from mechanical relaxation, and𝒗f\n𝑖 represents\nnormal motion due to tissue formation.\nThevelocity 𝒗m\n𝑖 isconstructedundertheassumptionofoverdampedmotion,appropriatefor\nslow cellular dynamics in viscous environments (Purcell,1977). In this regime, inertial effects\nareneglectedandvelocityisproportionaltothenetforceactingontheboundary 𝒓𝑖 (𝑡) suchthat\n𝒗m\n𝑖 = 1\n𝜂 𝑭(net)\n𝑖 ,(2)\n4\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigure2 –Overviewofthemathematicalmodel.(a)Tissuegrowthregion R andtime-dependent\ntissue area 𝛺(𝑡) showing differentiated cells (red) and the discrete cell-based tissue interface\n𝑆(𝑡). Tissue formation displaces interface cells perpendicularly in the direction ¯𝒏𝑖 with a\nvelocityproportionaltothelocalcelldensity 𝑞𝑖. (b)Mechanicalforcediagramatcellboundary\n𝒓𝑖 showing two restoring forces from springs on either side and the normal reaction force as\nwell as the net force. (c-f) Examples of discrete cells undergoing proliferation, asymmetric\nproliferation, apoptosis, and differentiation over a time stepΔ𝑡. Highlighted cells undergo the\nindicated process, and red cells denote differentiated cells as shown in (a).\nwhere 𝜂 is the viscous drag coefficient. The net force is composed of three contributions\n(Figure 2b): two restoring forces,𝑭(+)\n𝑖 and 𝑭(−)\n𝑖 , from the cells on either side, and a reaction\nforce 𝑭(n)\n𝑖 from the existing substrate. The restoring forces act along the cell bodies, with\nmagnitudesdeterminedbyaprescribedforcelaw 𝑓 (ℓ𝑖 (𝑡)) whichmaybeHookeanornonlinear,\nand ℓ𝑖 (𝑡) = ∥ 𝒓𝑖 (𝑡) − 𝒓𝑖−1(𝑡) ∥ is the cell length (Murray et al.,2012; Murphy et al.,2019).\nhe reaction force is defined to oppose the normal component of the restoring forces, thereby\nrestricting mechanical relaxation to the tangent direction of the interfaceˆ𝝉𝑖. Consequently, the\nnet force is given by the tangential projection\n𝑭(net)\n𝑖 =\nh\u0010\n𝑭(+)\n𝑖 + 𝑭(−)\n𝑖\n\u0011\n· ˆ𝝉𝑖\ni\nˆ𝝉𝑖. (3)\nThevelocitycontribution 𝒗f\n𝑖,modelstheperpendiculardisplacementofcellsinthedirection\n¯𝒏𝑖 as they deposit ECM at the interface (Figure2a). Compatibility considerations requires that\n5\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nthe magnitude of this velocity contribution to satisfy\n∥𝒗f\n𝑖 ∥=𝑘 f 𝑞𝑖,(4)\nwhere𝑞 𝑖 =1/ℓ 𝑖 is the density of the𝑖th cell (Buenzli,2015; Alias and Buenzli,2017; Hegarty-\nCremeretal., 2021;Kubaetal., 2026). Inthisdiscreteframework,perpendiculardisplacements\nof the cell segments can cause them to overlap or separate, which leads to a disconnected\ninterface. To combat this issue we implement a predictor-corrector algorithm to calculate𝒗f\n𝑖\ndepending on whether neighboring segments intersect or separate after being displaced. A full\ndescription of this algorithm is described in detail inKuba et al.(2026).\nInadditiontothemechanicalandformationprocesses,themodelincorporatesfourstochastic\ncell events:\n• Proliferation(Symmetricdivision): acelldividesintotwoactivetissue-formingdaugh-\nter cells that both lie on the tissue interface (Figure2c).\n• Proliferation (Asymmetric division): a cell divides into two daughter cells where one\nremains an active tissue-forming cell and the other differentiates and embeds into the\nECM (Figure2d).\n• Apoptosis: a cell dies and is then removed from the tissue interface (Figure2e).\n• Differentiation: a cell differentiates and embeds into the ECM (Figure2f).\nA key novelty of the framework is that it explicitly records the positions of differentiated\ncells (marked in red in Figure2a,d,f) during asymmetric division and differentiation, enabling\nsingle-cell resolution comparisons between model predictions and experimental observations.\nAlthough asymmetric division and direct differentiation both generate differentiated cells, they\nrepresent distinct cellular mechanisms. While these processes are not easily distinguishable in\nstatic experimental observations, they can be independently controlled and analysed within the\ncomputational framework.\nThese events are assumed to occur over a sufficiently small time intervalΔ𝑡, such that at\nmost one event takes place within[𝑡, 𝑡 + Δ𝑡) (Baker et al.,2019; Murphy et al.,2020). Each\ncell may undergo proliferation (either symmetric division or asymmetric division), apoptosis,\nor differentiation with constant rates (probabilities per unit time)𝑃sym, 𝑃asym, 𝐴, and 𝐷, re-\nspectively. Stochastic implementation proceeds by generating three random numbers from a\nuniformdistributiontodetermine: (i)whetheraneventoccurs, (ii)whicheventoccurs, and(iii)\nwhich cell undergoes the event (Gelman et al.,2013; Baker et al.,2019; Murphy et al.,2020).\nImplementation details are provided inA.1 and open-source code is available atKuba et al.\n(2026).\n6\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\n2.2 Connecting model parameters to experimental observations\nTo enable comparison between model-generated and experimental tissue compositions, model\nparameters were calibrated to reproduce experimentally measured differentiated cell densities.\nIn the model, differentiated cells arise at the moving tissue interface through either direct\ndifferentiation (𝐷) or asymmetric division (𝑃asym), such that the total differentiation rate is\n𝐷tot = 𝐷 + 𝑃asym. (5)\nFollowingBuenzli(2015), experimentally consistent differentiated cell densities𝜌 are obtained\nby approximating\n𝐷tot = 𝑘f 𝜌, (6)\nwhere 𝑘f is the tissue formation rate. For each case study, experimentally estimated values of𝜌\nand 𝑘f (Section 2.3) were substituted into Eq. (6) to determine𝐷tot.\nTo characterize the relative contribution of differentiation pathways, we define\n𝛼 = 𝐷\n𝐷tot\n,(7)\nrepresenting the fraction of differentiation occurring through direct differentiation. Since direct\ndifferentiationremovesanactivetissue-formingcellfromtheinterface,increasing𝛼reducesthe\ndensity of active cells and influences interface morphology via the normal velocity formulation\n(Eq. (4)).\nThe evolution of the active tissue-forming cell population𝑁is described by\nd𝑁\nd𝑡 =𝛽𝑁,(8)\nwhere𝛽=𝑃 sym −𝐴−𝐷 is the net growth rate. Here,𝛽 > 0 corresponds to a proliferative\npopulation and𝛽 < 0 to a declining population. Together,𝛼 and 𝛽 provide independent control\nover differentiation mode and population dynamics in the simulations\n2.3 Experimental data\nPreviouslypublishedexperimentalmeasurementsofdifferentiatedcelldensityandextracellular\nmatrix (ECM) formation rates were used to calibrate model parameters and ensure consistency\nbetween simulated and observed tissue compositions. Two experimental contexts were con-\nsidered: (i) tissue engineering cultures within 3D-printed constructs, and (ii) bone formation\nwithin cortical pores.\nFor tissue engineering applications, data were drawn from in vitro studies of MC3T3-E1\nosteoblasticcellsculturedinthin3D-printedporousscaffolds(Buenzlietal., 2020;Lanaroetal.,\n2021; Browning et al.,2021). In those experiments, cells were seeded within square pores of\n7\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nside length200–600μm and cultured for up to 28 days. Tissue growth was monitored using\nDAPI staining of cell nuclei and actin filaments, with ethidium homodimer staining indicating\ntheabsenceofcelldeath(Buenzlietal., 2020). Imageswereacquiredatmultipletimepoints(4,\n7, 10, 14, 18, and 28 days), from which differentiated cell density, cell orientation, tissue area\ncoverage, and pore infill duration were quantified (Lanaro et al.,2021; Browning et al.,2021).\nAlthoughorientation,coverage,andinfilltimevariedwithporesize,parameterestimationacross\nallporesindicatedthatatlatetime( 𝑡 > 14days)thedensityofdifferentiatedcellsremainaround\n𝜌 ≈ 0.003 cells/μm2 (Browning et al.,2021). The corresponding ECM formation rate was\nestimated as𝑘f ≈ 87.84 μm2/day (Kuba et al.,2026).\nFor cortical bone formation, matrix apposition rate (MAR) and osteocyte density (𝜌) were\nobtainedfrompreviouslypublisheddouble-labelingexperimentsandimagingstudiesinanimal\nmodels and human bone (Metz et al.,2003; Hannah et al.,2010; Power et al.,2012). Reported\nMAR values in human cortical bone range between0.5–3μm/day (Jones,1974; Martin et al.,\n1998; Buenzli et al.,2014; Alias and Buenzli,2018). To ensure consistency with this range\n(target ≈ 1.5 μm/day), the model parameter𝑘f was set to30 μm2/day, assuming a typical\nosteoblast size of∼ 20 μm and using Eq. (4). Consistent with reported estimates (Alias and\nBuenzli, 2018; Parfitt, 1983; Hannah et al.,2010; Power et al.,2012), osteocyte density was\ntaken as𝜌 = 0.000625 cells/μm2 throughout the pore.\n3 Results\nWe apply the computational model to two case studies: (i) tissue growth within 3D-printed\nscaffolds (Section3.1) and (ii) bone formation within cortical pores (Section3.2). In the tissue\nengineering study, we investigate how differentiation mechanisms, cell population dynamics,\nand mechanical properties of cells shape tissue growth patterns, interface morphology, and\ndifferentiated cell orientations. In the bone formation case study we examine how stochasticity\ninfluences the distribution of differentiated cells and the structural symmetry of the pore.\n3.1 Tissue Engineering\nTo reflect the experimental setup described in Section2.3, the tissue interface was initialized\nin a 300 μm × 300 μm square with 𝑁 = 60 uniformly distributed cells of length20 μm,\napproximatingosteoblasticsize. FromEq.(6),weestimated 𝐷tot ≈ 0.26,andduetotheabsence\nof cell death we set𝐴 = 0. Since the fraction of differentiated cells (𝛼) that emerge from\ndirectdifferentiationandthenetgrowthrateofthetissue-formingpopulation( 𝛽)aredifficultto\nmeasure experimentally, we investigated their effects by varying these parameters.\nFigure 3a shows representative model-generated tissue compositions, with differentiated\ncells shown in red and the tissue interface indicated at the initial time (dark blue) and at 90%\npore fill (blue). Within the newly formed tissue, differentiated cells vary in orientation and size\n8\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigure 3– Example realisations of the tissue growth model for a series of control parameters\n𝛼and𝛽values. Simulations are initialised with a square geometry to emulate the cell culture\nexperiments with𝐷 tot =0. 003. (a) Realisations of the model for different values of𝛼 with\n𝛽 = 0. (b)Exampleimagesofexperimentsculturedwithing 300 μm × 300 μm,reproducedwith\npermissions fromBuenzli et al.(2020). (c) Solid lines are the average normalized area of void\nspace over time for different values of𝛽 with 𝛼 = 0.2, and black points are averages estimated\nfromexperimentsculturedinscaffoldsofthesamesize. Thesimulatedaveragesweregenerated\nfrom from a sample size of 50 realizations for each value of𝛽.\n9\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigure4–Simulationsofthestochasticmodelshowingtheeffectofcellstiffnessonembedded\ncell orientations. Each column shows a single realization for (a) low (𝑘=5 N/μm) (b)\nintermediate (𝑘 = 100 N/ μm) and (c) high (𝑘 = 500 N/ μm) stiffness. The top row displays\ntissue growth with differentiated cells (red) and the interface at𝑡 = 0, 5, and 10 days, using a\ncolor gradient from dark blue (early) to light blue (late). Subsequent rows show the frequency\nof cell orientations measured relative to the x-axis and binned into 24 intervals over[0◦ , 360◦ ).\n10\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigure 5– Comparison of experimental and simulated cell orientation distributions. (a) Mea-\nsuredcellorientationsfromcellcultureexperimentsondays 4, 7, and 10ina 300 μm × 300 μm\n3D-printed square porous scaffold reported byLanaro et al. (2021) (reproduced with per-\nmissions). (b-d) Cumulative cell orientation frequencies from 100 realizations for low\n(𝑘 = 5 N/ μm), intermediate (𝑘 = 100 N/ μm) and high (𝑘 = 500 N/ μm) cell stiffness at\nsimulation times that correspond to the experimental data.\nbutremainapproximatelyuniformlydistributed,despitetheincreasingdensityoftissue-forming\ncells as the interface evolves. To investigate the influence of differentiation modes, in Figure3a\nwevaried 𝛼andset 𝛽 = 0sothat 𝑃sym, 𝑃asym,and 𝐷 balancetomaintainaconsistentpopulation\nof tissue-forming cells. Comparison with experimental observations (Figure3b) indicates that\n𝛼 = 0.2 provides the closest agreement between model predictions and experiments, whereas\nlarger 𝛼 values yield a rougher and less consistent interface morphology. Figure3c shows the\neffectofadynamictissue-formingcellpopulationonthetimerequiredtoachieveporeinfilling.\nIn these simulations,𝛼 was fixed at0.2, and growing (𝛽 > 0), constant (𝛽 = 0) and shrinking\n(𝛽 < 0) populations were considered. As in Figure3a, simulations were performed until 90%\npore fill, and the normalized tissue area, averaged from 50 realizations, was plotted over time\nfor each𝛽 value alongside the averaged experimental data for scaffolds of this size.\nFigure4examinestheinfluenceofcellmechanicalstiffnessondifferentiatedcellorientation.\nStiffnesslevelswerecategorizedaslow( 𝑘 = 5 N/μm),intermediate( 𝑘 = 100 N/μm),andhigh\n(𝑘 = 500 N/μm). Representativerealizationsandcorrespondingroseplotsshowthatorientation\ndistributions are initially anisotropic across all stiffness levels (𝑡 = 4 days), but anisotropy de-\n11\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\ncreasesovertimeforlowerstiffnessvalues. Comparisonofcumulativeorientationdistributions\nwith experimental measurements (Lanaro et al., 2021) indicates improved agreement at lower\nstiffness, although the associated interface morphologies do not reproduce the experimentally\nobserved rounding (Figure 5).\n3.2 Bone formation\nTo model bone tissue formation within irregular resorption cavities, the tissue interface was\ninitialized as a deterministically perturbed circle. Specifically,𝑁values were sampled from a\nuniform distribution and smoothed with respect to the angular coordinate𝜃 ∈ [ 0, 2𝜋) using the\nloessfunction to generate spatially correlated perturbations (Cleveland and Grosse,1991). The\nresulting smoothed values𝑢𝑖 ∈ [ 0, 1] were scaled by a perturbation amplitude˜𝑅 and added to a\ncircle of radius𝑅0 such that, at𝑡 = 0,\n𝑅𝑖 = 𝑅0 + 𝑢𝑖 ˜𝑅, 𝑖 = 1, . . . , 𝑁 . (9)\nExperimental evidence indicates that bone formation initiates once osteoblast density reaches a\nthreshold, at which point proliferation ceases (Lassen et al.,2017). Accordingly, proliferation\nrates were set to zero (𝑃sym = 𝑃asym = 0), so differentiation occurred solely through direct\ndifferentiation, yielding𝐷 = 𝑘f𝜌 ≈ 0.019. The apoptosis rate𝐴 was chosen to be lower than𝐷\nto prevent rapid cell loss.\nFigure 6b shows a representative realization for an artificial resorption cavity, with osteo-\ncytes marked in black and growth regions shaded according to formation time, consistent with\nexperimental images (Figure6a). Osteocyte density was quantified within regions formed over\na time interval of𝑇 = 7 days, corresponding to lamellar deposition in cortical bone. Summary\nstatistics (Figure6c) demonstrate that local osteocyte density remains consistent and closely\nmatches the prescribed density𝜌, despite variation in region shape and size over time.\nAlthoughtherealizationinFigure 6bshowssymmetricporeclosure,corticalporesmayalso\ndevelopasymmetricstructures. Asymmetrywasquantifiedusingthewallthicknessasymmetry\nratio W.Thasym, defined using the pair of opposing wall thicknesses with the greatest disparity\n(Figure6a), where wall thickness is measured as the distance between the inner and outer pore\nboundaries (Hegarty-Cremer et al.,2024):\nW.Thasym = W.Thmin\nW.Thmax\n.(10)\nFigure 6d shows representative realizations classified by W.Thasym, and Figure 6e summarizes\nosteocyte densities in the wall regions used to compute this ratio. These results indicate that\nosteocyte density is higher near thinner walls, and the most frequently observed asymmetries\nfall within the range W.Thasym ∈ [0.6,0.7).\n12\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigur\ne 6– Model simulations of bone formation in an artificial osteon. (a) Experimental image illustrating the measurement of W.Thmax and W.Thmin\nused to compute the W.Thasym ratio. (b) Representative model realization in an artificial osteon, with osteocytes marked in black. Pseudo-lamellae\nformed over a time interval𝑇 = 7 days are shaded according to formation time. (c) Box plots of embedded cell densities within wall regions formed\nover 𝑇 days from 500 realizations. Shading corresponds to the wall regions illustrated in (a). (d) Example model realizations spanning a range of\nW.Thasym values. (e) Summary statistics of embedded cell densities measured in rectangular regions nearW.Thmax (blue) andW.Thmin (red), grouped\nby asymmetry ratio bins. Example sampling regions are highlighted in the top-left panel of (d). In all model realizations death rate of cells is𝐴 = 0.01\nand total differentiation rate is𝐷tot = 𝑘f 𝜌 where 𝜌 = 0.000625 cells/μm2.\n13\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\n4 Discussion\nExperimental evidence indicates that biological tissue growth is strongly regulated by substrate\ngeometry and mechanobiological cues (Nelson et al.,2005; Rumpler et al.,2008; Bidan et al.,\n2012;LadouxandMége, 2017;Callensetal., 2020). Inconfinedgeometries,spatialconstraints\ncausetissue-formingcellstocrowdorspreadduringgrowth,generatingmechanicalstressesthat\ninturninfluencecellularprocessessuchasproliferationanddifferentiation(Nelsonetal., 2005;\nWerner et al.,2017; Gudipaty et al.,2017). However, disentangling the mechanistic interplay\nbetweengeometry,mechanics,andcellularprocessesintissuegrowthremainschallengingwhen\nrelyingsolelyonexperimentaltissuesamples,eventhoughsuchsamplescontainfingerprintsof\ntheprocessesthatshapedthem. Toaddressthislimitation,weproposeacomputationalmodelof\ntissuegrowththatgeneratessynthetictissueswithsingle-cellresolution,enablingdirectcontrol\nof governing mechanisms and systematic comparison with experimental data to isolate their\nindividual contributions.\nDuringtissuegrowth,cellsmayundergoproliferation,differentiation,andapoptosis,withthe\nrelative contributions of these processes varying across tissue types. Histological staining can\nreveal the presence or absence of certain processes, but it cannot resolve the specific processes\nthat occurred during growth. In tissue engineering experiments (Figure1a), DAPI staining\nindicates that there is a substantial amount of cell proliferation during growth, while ethidium\nhomodimer staining indicates the absence of cell death. Model simulations indicate that, given\nthe incorporated mechanisms, closer agreement with experimental observations and smoother,\nmore consistent tissue interfaces are achieved when the majority of differentiated cells arise\nfrom asymmetric division (𝛼 = 0.2) rather than direct differentiation (Figure3a). In addition,\ncomparisonofnormalizedtissueareaovertimeindicatesthatanincreasingpopulationoftissue-\nforming cells more accurately captures late-stage growth dynamics (Figure3c), consistent with\nprevious quantitative and modeling studies (Buenzli et al.,2020; Browning et al.,2021).\nInthesetissueengineeredexperiments,differentiatedcellsareassumedtobecomequiescent\nonceembeddedwithintheECM.Consequently,differentiatedcellorientationsprovidearecord\nof the position and shape of the tissue interface at the time of differentiation. Previous mod-\neling studies of evolving tissue interfaces have demonstrated a relationship between interface\nmorphologyandthestrengthofmechanicalinteractionsbetweencells(AliasandBuenzli, 2017;\nKubaetal., 2026). Consistentwiththis, modelsimulationsshowthatintermediatecellstiffness\nproduces tissue interfaces that more closely resemble the experimentally observed rounding\n(Figure 4), whereas low stiffness yields orientation distributions with reduced anisotropy that\nbetter match experimental measurements (Lanaro et al.,2021) (Figure5). This discrepancy\nsuggests that additional mechanisms may contribute to tissue growth mechanics beyond those\ncaptured by cell–cell mechanical interactions alone.\nA plausible contributor to these additional mechanisms is surface tension. In their experi-\nmental study,Lanaro et al.(2021) observed groups of cells spanning perpendicular fibers and\n14\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nFigure 7– Example realisations of the tissue growth model in different initial geometries. The\ngeometries are motivated by previously conducted cell culture experiments ofRumpler et al.\n(2008) (triangle),Bidan et al.(2013) (cross), andVandenHeuvel et al.(2023) (wavy).\nproposed that cells actively seek nearby attachment points to maintain surface tension during\ngrowth. This behavior is consistent with purse–string–based closure mechanisms observed in\nwound healing (Ravasio et al.,2015) and with evidence that surface tension plays a central\nroleintissueorganizationandmechanicalregulation(Bidanetal., 2016). Inthecurrentmodel,\nsurfacetension–relatedforcesactingperpendiculartotheinterfacearenotrepresentedduetothe\ndefinitionofthereactionforce 𝑭(n)\n𝑖 . Thesefindingsunderscoretheimportanceofincorporating\nsurface tension–mediated effects into tissue growth models, as neglecting them may obscure\nkeymechanicalcontributionstointerfaceevolution. Althoughthiscasestudyfocusedonsquare\ngeometries,thesameanalysiscanbeextendedtootherporeshapes(Figure 7)tofurtherexamine\nthe role of surface tension in tissue growth.\nInSection 3.2, weinvestigatedhowstochasticityinfluencesthespatialdistributionofdiffer-\nentiated cells and the resulting growth patterns during bone formation within cortical pores. In\nthis setting, a finite population of osteoblasts deposits new matrix at the pore surface without\nproliferation, differentiating into osteocytes that become embedded within the mineralized tis-\nsue. Experimental studies report that osteocyte density decreases with distance from the initial\npore boundary (Power et al.,2012), although differentiation rates cannot be directly measured.\nConsistent with previous modeling work (Buenzli,2015), our simulations show that osteocyte\ndensity within newly formed lamellar regions remains consistent with the prescribed constant\n15\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nvalue𝜌, despite variations in regional shape and size (Figure 6b,c). This suggests that dy-\nnamic differentiation rates could be inferred from single experimental snapshots by quantifying\nosteocyte densities within lamellae and rearranging Eq. (6), providing a basis for incorporat-\ning curvature-, stress-, or property-dependent differentiation rules (e.g.,𝐷 (𝜅, 𝜎, . . .)) into the\nmodeling framework.\nCorticalboneformationoccursaroundHaversiancanalswhosepositionsareoftenasymmet-\nric(Robling, 1998;SchnitzlerandMesquita ,2013;Cookeetal., 2022),motivatinginvestigation\ninto the mechanisms underlying such asymmetries. Previous modeling studies have attributed\nasymmetricgrowthtobiologicalheterogeneities, suchasspatialvariationsinmatrixdeposition\nrates (Hegarty-Cremer et al.,2024). In contrast, our analysis (Figure6d,e) demonstrates that\nasymmetriescanarisepurelyfromstochasticdifferentiationandapoptosisofbone-formingcells.\nThediscreteformulationofthemodelcapturesdensityfluctuationsassociatedwiththeseevents,\nleading to locally reduced normal velocities (Eq. (4)) and emergent interface irregularities.\nThese findings suggest that at least part of the variability observed in cortical pore morphology\nmayresultfromintrinsic stochasticcellularprocesses, althoughfurtherexperimentalvalidation\nis required.\n5 Conclusions and future work\nExperimental evidence demonstrates that tissue growth is governed by the interplay between\ngeometry,mechanics,andbiologicalprocesses(Nelsonetal., 2005;Rumpleretal., 2008;Bidan\net al.,2012; Ladoux and Mége,2017; Callens et al.,2020). Disentangling these mechanisms\nexperimentally remains challenging, as mechanical properties are difficult to measure (Tee\net al.,2009) and tissue samples provide only static snapshots of a dynamic process. However,\nthe material deposited behind the advancing tissue interface retains structural signatures of\nthe processes that shaped tissue architecture. Here, we introduced a computational model\nthatsimulatestissuegrowthatsingle-cellresolutionandgeneratessynthetictissuecompositions\ndirectlycomparablewithexperimentalobservations. Throughtwocasestudies,wedemonstrated\nthattheframeworkcan(i)identifythecontributionofspecificmechanismstotissuemorphology,\n(ii)estimatedynamicratesfromsingle-time-pointdata,and(iii)generatetestablehypothesesfor\nthe origins of stochastic growth patterns. These results highlight the utility of model-generated\ntissue compositions as a complementary tool for investigating tissue growth dynamics and\nenabling in silico hypothesis testing.\nFutureextensionscouldfocusontailoringtheframeworktospecifictissuetypesbyincorpo-\nratingbiologicallyinformedprocessrules. Forexample,inboneformation,differentiationrates\ncould be extended to depend on interface curvature or mechanical stress and calibrated using\nexperimental estimates. Similar extensions could be applied to proliferation, apoptosis, and\nthe tissue formation rate𝑘f, consistent with previous studies in which𝑘f depended on interface\nradiusorporosity(AliasandBuenzli, 2018), andcelleventratesdependedoncelllength(Mur-\n16\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nphy et al.,2020). Parameter estimation could then be performed by comparing synthetic and\nexperimental datasets using statistical inference techniques. Additionally, the model’s ability\nto track embedded cell positions enables spatial statistical analyses, such as point process or\npair-correlationmethods,providingaquantitativeframeworkfortestingcompetingmechanistic\nhypotheses and guiding future experimental design.\nAcknowledgements\nThisresearchwassupportedbytheAustralianResearchCouncil(DP190102545,DP230100025).\nSK and PRB acknowledge support from the Max Planck Queensland Centre for the Materials\nScienceofExtracellularMatrices. WethankBrennaDevlinandMariaA.Woodruffforproviding\nthe experimental image presented in Figure 1a.\nAppendix A Mathematical Model\nA.1 Implementation of stochastic cellular processes in mathematical\nmodel\nTo simplify notation, we write the probability that the𝑖th cell proliferates is𝑃sym 𝑖 (symmetric\ndivision). Similarly for asymmetric proliferation, apoptosis, and differentiation, we write the\nprobabilities as 𝑃asym 𝑖, 𝐴𝑖, and 𝐷𝑖. The algorithm at each numerical time stepΔ𝑡 has been\nimplemented as follows:\n1. Generate 𝑝1, 𝑝2, 𝑝3 from a uniform distribution in[0, 1].\n2. If 𝑝1 ≤ Í𝑁\n𝑖=1(𝑃sym 𝑖 + 𝑃asym 𝑖 + 𝐴𝑖 + 𝐷𝑖)Δ𝑡, where𝑁 isthetotalnumberofcellsthatform\nthe interface, then one of the cells is experiencing one of the cellular processes over the\ntime interval[𝑡 + Δ𝑡). Otherwise continue to the next simulation time step.\n3. Using 𝑝2, we determine which cellular process occurred under the following conditions:\n• If 𝑝2 ≤\nÍ𝑁\n𝑖=1 𝑃sym 𝑖\nÍ𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) then the occurring cellular process is proliferation.\n•Else if\nÍ𝑁\n𝑖=1 𝑃sym𝑖\nÍ𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) < 𝑝 2 ≤\nÍ𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 )Í𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) then the occurring\ncellular process is asymmetric proliferation.\n•Else if\nÍ𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 )Í𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) < 𝑝 2 ≤\nÍ𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖)Í𝑁\n𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) then the occurring\ncellular process is apoptosis.\n•Else the occurring cellular process is differentiation.\n4. Using𝑝 3,wedeterminewhichcellexperiencedthecellularprocess. Forexample,assum-\ningthatthepreviousstepresultedinproliferation(symmetricdivision),weassertthatthe\n17\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\n𝑗th cell proliferates when the conditions below are met,\n0< 𝑝 3 ≤\nÍ 𝑗\n𝑖=1 𝑃sym𝑖\nÍ𝑁\n𝑖=1 𝑃sym𝑖\n, 𝑗=1,\nÍ 𝑗−1\n𝑖=1 𝑃sym𝑖\nÍ𝑁\n𝑖=1 𝑃sym𝑖\n< 𝑝 3 ≤\nÍ 𝑗\n𝑖=1 𝑃sym𝑖\nÍ𝑁\n𝑖=1 𝑃sym𝑖\n,2≤ 𝑗 ≤ 𝑁 .\nA similar process is followed in the case that the previous step results in apoptosis,\nembedment, and simultaneous proliferation and embedment.\n5. Based on cell event, we either add or remove cell boundaries from the interface, linearly\ninterpolating the boundary nodes over the region where the𝑗th cell spanned as described\nin Section2.1.\nReferences\nPoujade, M., Grasland-Mongrain, E., Hertzog, A., Silberzan, P.,Collective migration of an\nepithelialmonolayerinresponsetoamodelwound.,CellBiology,2007;104: 15988–15993.\ndoi: 10.1073/pnas.0705062104.\nRolli,C.G.,Nakayama,H.,Yamaguchi,K.,Spatz,J.P.,Kemkemer,R.,Nakanishi,J., Switchable\nadhesive substrates: Revealing geometry dependence in collective cell behavior., Biomateri-\nals, 2012; 33: 2409–2418. doi:10.1016/j.biomaterials.2011.12.012.\nMartin, R.B., Burr, D.B., Sharkey, N.A., Fyhrie, D.P.Skeletal Tissue Mechanics., New York:\nspringer, 1998.\nLowengrub,J.S.,Frieboes,H.B.,Jin,F.,Chuang,Y-L.,Li,X.,Macklin,P.,Wise,S.M.,Cristini,\nV.,Nonlinearmodellingofcancer: bridgingthegapbetweencellsandtumours.,Nonlinearity,\n2009; 23: R1. doi:10.1088/0951-7715/23/1/R01 .\nWozniak,P.,Elhaj,A.J., TissueEngineeringUsingCeramicsandPolymers,Chapter14-Bone\nregeneration and repair using tissue engineering., Woodhead Publishing, 2007.\nHollister, S.J.,Porous scaffold design for tissue engineering, Nature Materials, 2005; 4: 518–\n524. doi:10.1038/nmat1421.\nO’Brien, F.J.,Biomaterials & scaffolds for tissue engineering, Materials Today, 2011; 14:\n88–95. doi:10.1016/S1369-7021(11)70058-X.\nBidan,C.M.,Kommareddy,K.P.,Rumpler,M.,Kollmannsberger,P.,Fratzl,P.,Dunlop,J.W.C.,\nGeometryasaFactorforTissueGrowth.,AdvancedHealthcareMaterials,2013;2: 186–194.\ndoi: 10.1002/adhm.201200159.\n18\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nRumpler, M., Woesz, A., Dunlop, J.W.C, van Dongen, J.T, Fratzl, P.,The effect of geometry on\nthree-dimensional tissue growth., Journal of The Royal Society Interface, 2008; 27: 1173–\n1180. doi: 10.1098/rsif.2008.0064.\nBidan, C., Kommareddy, K., Rumpler, M., Kollmannsberger, P., Bréchet, Y., Fratzl, P., Dunlop,\nJ.,How Linear Tension Converts to Curvature: Geometric Control of Bone Tissue Growth.,\nPloS one, 2012; 7: e36336. doi: 10.1371/journal.pone.0036336.\nSchamberger, B., Ziege, R., Anselme, K., Ben Amar, M., Bykowski, M., Castro, A.P.G.,\nCipitria, A., Coles, R.A., Dimova, R., Eder, M., Ehrig, S., Escudero, L.M., Evans, M.E.,\nFernandes, P.R., Fratzl, P., Geris, L., Gierlinger, N., Hannezo, E., Iglič, A., Kirkensgaard,\nJ.J.K.,Kollmannsberger,P.,Kowalewska,Ł.,Kurniawan,N.A.,Papantoniou,I.,Pieuchot,L.,\nPires, E.F., Roschger, A., Bidan, C.M., Dunlop, J.W.C.,Curvature in Biological Systems: Its\nQuantification, Emergence, and Implications across the Scales., Advanced Materials, 2023;\n35: 2206110. doi:10.1002/adma.202206110.\nCallens,S.J.P.,Uyttendaele,R.J.C.,Fratila-Apachitei,L.E.,Zadpoor,A.A., Substratecurvature\nas a cue to guide spatiotemporal cell and tissue organization., Biomaterials, 2020; 232:\n119739. doi:10.1016/j.biomaterials.2019.119739.\nMcNamara,L.E.,McMurray,R.J.,Dalby,M.J., NanotopographicalControlofStemCellDiffer-\nentiation., Journal of Tissue Engineering, 2010; 1: 120623. doi:10.4061/2010/120623.\nDobbenga, S., Fratila-Apachitei, L.E., Zadpoor, A.A.,Nanopattern-induced osteogenic differ-\nentiation of stem cells – A systematic review., Acta Biomaterialia, 2016; 46: 3–14. doi:\n10.1016/j.actbio.2016.09.031.\nNelson, C.M., Jean, R.P., Tan, J.L., Liu, W.F., Sniadecki, N.J., Spector, A.A., Chen, C.S.,\nEmergentpatternsofgrowthcontrolledbymulticellularformandmechanics.,Proceedingsof\ntheNationalAcademyofSciences,2005;102: 11594–11599.doi: 10.1073/pnas.0502575102.\nDischer, D.E., Janmey, P., Wang, Y-L.,Tissue Cells Feel and Respond to the Stiffness of Their\nSubstrate., American Association for the Advancement of Science, 2005; 310: 1139–1143.\ndoi: 10.1126/science.1116995.\nLadoux,B.,Mége,R.M., Mechanobiologyofcollectivecellbehaviours.,NatureReviewsMolec-\nular Cell Biology, 2017; 18: 743–757. doi:10.1038/nrm.2017.98.\nLim, C.T., Bershadsky, A., Sheetz M.P.,Mechanobiology., Journal of The Royal Society Inter-\nface, 2010; 7: S291–S293. doi:10.1098/rsif.2010.0150.focus.\nAlberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K., Walter, P.,Molecular Biology of the\nCell, 4th edition., New York: Garland Science, 2002.\n19\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nKurniawan, N.A., Bouten, C.V.C.,Mechanobiology of the cell–matrix interplay: Catching a\nglimpseofcomplexityviaminimalisticmodels,ExtremeMechanicsLetters,2018;20: 59–64.\ndoi: 10.1016/j.eml.2018.01.004.\nSkedros,J.G.,Holmes,J.L.,Vajda,E.G.,Bloebaum,R.D., Cementlinesofsecondaryosteonsin\nhumanbonearenotmineral-deficient: Newdatainahistoricalperspective.,TheAnatomical\nRecord, 2005; 286: 781–803. doi:10.1002/ar.a.20214.\nRoschger,P.,Paschalis,E.P.,Fratzl,P.,Klaushofer,K., Bonemineralizationdensitydistribution\nin health and disease., Bone, 2008; 42: 456–466. doi:10.1016/j.bone.2007.10.021.\nVatsa, A., Breuls, R.G., Semiens, C.M., Salmon, P.L., Smit, T.H., Klein-Nulend, J.,Osteocyte\nmorphology in fibula and calvaria — Is there a role for mechanosensing?, Bone, 2008; 43:\n452–458. doi:10.1016/j.bone.2008.01.030.\nAscenzi, M., Gills, J., Lomovstev, A., Orientation of collagen at the osteocyte lacunae\nin human secondary osteons., Journal of Biomechanics, 2008; 41: 3426–3435. doi:\n10.1016/j.jbiomech.2008.09.010.\nvan Hove, R.P., Nolte, P.A., Vatsa, A., Semiens, C.M., Salmon, P.L., Smit, T.H., Klein-Nulend,\nJ., Osteocyte morphology in human tibiae of different bone pathologies with different bone\nmineral density — Is there a role for mechanosensing?, Bone, 2009; 45: 321–329. doi:\n10.1016/j.bone.2009.04.238.\nBritz, H.M., Carter, Y., Jokihaara, J., Leppanen, O.V., Jarvinen, T.L.N., Belev, G., Cooper,\nD.M.L.,Prolonged unloading in growing rats reduces corticalosteocyte lacunar density and\nvolume in the distal tibia., Bone, 2012; 51: 913–919. doi:10.1016/j.bone.2012.08.112.\nMader, K.S., Schneider, Philipp, S., Muller, R., Stampanoni, M.,A quantitative framework for\nthe 3D characterization of the osteocyte lacunar system., Bone, 2013; 57: 142–154. doi:\nhttps://doi.org/10.1016/j.bone.2013.06.026.\nLanaro, M., Mclaughlin, M.P., Simpson, M.J., Buenzli, P.R., Wong, C.S., Allenby,\nM.C., Woodruff, M.A., A quantitative analysis of cell bridging kinetics on a scaf-\nfold using computer vision algorithms., Acta Biomaterialia, 2021; 136: 429–440. doi:\nhttps://doi.org/10.1016/j.actbio.2021.09.042.\nRuffoni, D., Fratzl, P., Roschger, P., Klaushofer, K., Weinkamer, R.,The bone mineralization\ndensity distribution as a fingerprint of the mineralization process., Bone, 2007; 40: 1308–\n1319. doi:10.1016/j.bone.2007.01.012.\nBuenzli, P.R., Osteocytes as a record of bone formation dynamics., Journal of Theoretical\nBiology, 2015; 364: 418–427. doi:10.1016/j.jtbi.2014.09.028.\n20\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nCox, B.N., Purohit, P.K., White, S.N.Strain fields and solitary strain waves as determining\nfactors for the cross-sectional geometry of mouse incisor enamel., Journal of the Mechanics\nand Physics of Solids, 2015; 193: 105840. doi:10.1016/j.jmps.2024.105840.\nGuyot, Y., Papantoniou, I., Chai, Y.C., Van Bael, S., Schrooten, J., Geris, L.,A computational\nmodel for cell/ECM growth on 3D surfaces using the level set method: a bone tissue engi-\nneering case study., Biomechanics and Modeling in Mechanobiology, 2014; 13: 1361–1371.\ndoi: 10.1007/s10237-014-0577-5.\nCallens, S.J.P., Fan, D., van Hengel, I.A.J., Minneboo, M., Díaz-Payno, P.J., Stevens, M.M.,\nFratila-Apachitei, L.E., Zadpoor, A.A.,Emergent collective organization of bone cells in\ncomplex curvature fields., Nature Communications, 2023; 14: 855. doi:10.1038/s41467-\n023-36436-w.\nAmbrosi, D., Guana, F.,Stress-Modulated Growth., Mathematics and Mechanics of Solids,\n2007; 12: 319–342. doi:10.1177/1081286505059739.\nDunlop,J.W.C.,Fischer,F.D.,Gamsjäger,E.,Fratzl,P., Atheoreticalmodelfortissuegrowthin\nconfined geometries., Journal of the Mechanics and Physics of Solids, 2010; 58: 1073–1087.\ndoi: 10.1016/j.jmps.2010.04.008.\nGoriely, A.,The Mathematics and Mechanics of Biological Growth., Springer New York, NY,\n2017.\nJoly, P., Duna, G.N., Schone, M., Welzel, P.B., Fruedenberg, U., Werner, C., Petersen, A.,\nGeometry-Driven Cell Organization Determines Tissue Growths in Scaffold Pores: Conse-\nquences for Fibronectin Organization., PLOS One, 2013; 8: e73545. doi:10.1371/jour-\nnal.pone.0073545.\nGahffarizadeh,A.,Heiland,R.,Friedman,S.H.,Mumenthaler,S.M.,Macklin,P., PhysiCell: An\nopensourcephysics-basedcellsimulatorfor3-Dmulticellularsystems.,PLOSComputational\nBiology, 2018; 14: e1005991. doi:10.1371/journal.pcbi.1005991.\nAlias, M.A., Buenzli, P.R., Modeling the Effect of Curvature on the Collective Behav-\nior of Cells Growing New Tissue., Biophysical Journal, 2017; 112: 193–204. doi:\n10.1016/j.bpj.2016.11.3203.\nHegarty-Cremer,S.G.D.,Simpson,M.J.,Andersen,T.L.,Buenzli,P.R., Modellingcellguidance\nand curvature control in evolving biological tissues., Journal of Theoretical Biology, 2021;\n520: 110658. doi:10.1016/j.jtbi.2021.110658.\nKuba, S., Simpson, M.J., Buenzli, P.R., A mathematical model of curvature con-\ntrolled tissue growth incorporating mechanical cell interactions, 2026; BioRxiv. doi:\n10.64898/2026.03.10.710423.\n21\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nBaker, R.E., Parker, A., Simpson, M.J.,A free boundary model of epithelial dynamics., Journal\nof Theoretical Biology, 2019; 481: 61–74. doi:10.1016/j.jtbi.2018.12.025.\nMurphy, R.J., Buenzli, P.R., Baker, R.E., Simpson, M.J., A one-dimensional individual-\nbased mechanical model of cell movement in heterogeneous tissues and its coarse-\ngrained approximation., Proceedings of the Royal Society A, 2019; 475: 20180838. doi:\n10.1098/rspa.2018.0838.\nBuenzli, P.R., Kuba, S., Murphy, R.J., Simpson, M.J.,Mechanical cell interactions on curved\ninterfaces., Bulletin of Mathematical Biology, 2025; 87: 29. doi: 10.1007/s11538-024-\n01406-w.\nBrown,P.J.,Green,E.F.,Binder,B.J.,Osborne,J.M., Competingmechanismsforthebucklingof\nanepithelialmonolayeridentifiedusingmulticellularsimulation.,MathematicalBiosciences,\n2025; 380: 109367. doi:https://doi.org/10.1016/j.mbs.2024.109367.\nMurray, P.J., Edwards, C.M., Tindall, M.J., Maini, P.K.,From a discrete to a continuum\nmodel of cell dynamics in one dimension., Physical Review E, 2009; 80: 031912. doi:\n10.1103/PhysRevE.80.031912.\nMurray, P.J., Edwards, C.M., Tindall, M.J., Maini, P.K.,Classifying general nonlinear force\nlaws in cell-based models via the continuum limit., Physical Review E, 2012; 85: 021921.\ndoi: 10.1103/PhysRevE.85.021921.\nGelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B.,Bayesian Data\nAnalysis, Third Edition, Chapman and Hall/CRC, 2013.\nCleveland, W.S., Grosse, E.,Computational methods for local regression, Statistics and Com-\nputing, 1991; 1(1), 47-62. doi:10.1007/BF01890836.\nPurcell, E.J.,Life at low Reynolds number, American Journal of Physics 1977; 45: 3–11. doi:\n10.1142/9789814434973_0004.\nMurphy, R.J., Buenzli, P.R., Baker, R.E., Simpson, S.J., Mechanical Cell Competition in\nHeterogeneous Epithelial Tissues., Bulletin of Mathematical Biology, 2020; 82: 130. doi:\n10.1007/s11538-020-00807-x.\nBuenzli, P.R., Lanaro, M., Wong, C.S., McLaughlin, M.P., Allenby, M.C., Woodruff, M.A.,\nSimpson, M.J., Cell proliferation and migration explain pore bridging dynamics in 3D\nprinted scaffolds of different pore size., Acta Biomaterialia, 2020; 114: 285–295. doi:\n10.1016/j.actbio.2020.07.010.\nBrowning, A.P., Maclaren, O.J., Buenzli, P.R., Lanaro, M., Allenby, M.C., Woodruff,\nM.A., Simpson, M.J., Model-based data analysis of tissue growth in thin 3D\n22\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nprinted scaffolds., Journal of Theoretical Biology, 2021; 528: 110852. doi:\nhttps://doi.org/10.1016/j.jtbi.2021.110852.\nParfitt, A.M.,The physiological and clinical significance of bone histomorphometric data., In\nRecker R.R. (Ed.), Bone histomorphometry: Techniques and interpretation., 1983.\nAlias, M.A., Buenzli, P.R.,Osteoblasts infill irregular pores under curvature and porosity\ncontrols., Biomechanics and Modeling in Mechanobiology, 2018; 17: 1357–1371. doi:\n10.1007/s10237-018-1031-x.\nMetz, L.N., Martin, R.B., Turner, A.S.,Histomorphometric analysis of the effects of osteo-\ncyte density on osteonal morphology and remodeling., Bone, 2003; 33: 753–759. doi:\n10.1016/S8756-3282(03)00245-X.\nHannah, K.M., Thomas, C.D.L., Clement, J.G., Peele, A.G.,Bimodal distribution of osteocyte\nlacunarsizeinthehumanfemoralcortexasrevealedbymicro-CT.,Bone,2010;47: 866–871.\ndoi: https://doi.org/10.1016/j.bone.2010.07.025.\nPower, J., Doube, M., van Bezooijen, R.L., Loveridge, N., Reeve, J.,Osteocyte recruitment\ndeclines as the osteon fills in: Interacting effects of osteocytic sclerostin and previous hip\nfracture on the size of cortical canals in the femoral neck., Bone, 2012; 50: 1107–1114. doi:\nhttps://doi.org/10.1016/j.bone.2012.01.016.\nJones S.P.,Secretory territories and rate of matrix production of osteoblasts., Calcified Tissue\nResearch, 1974; 14: 309–315. doi:10.1007/BF02060305.\nMartin, R.B., Burr, D.B., Sharkey, N.A.,Skeletal Tissue Mechanics., Springer New York, NY,\n1998.\nBuenzli, P.R., Pivonka, P., Smith, D.W.,Bone refilling in cortical basic multicellular units:\ninsights into tetracycline double labelling from a computational model., Biomechanics and\nModeling in Mechanobiology, 2014; 13: 185–204. doi:10.1007/s10237-013-0495-y.\nBidan, C.M., Kollmannsberger, P., Gering, V., Ehrig, S., Joly, P., Petersen, A., Vogel, V., Fratzl,\nP., Dunlop, J.W.C.,Gradual conversion of cellular stress patterns into pre-stressed matrix\narchitecture during in vitro tissue growth., Journal of The Royal Society Interface, 2016; 13:\n20160136. doi:10.1098/rsif.2016.0136.\nVandenHeuvel,D.J.,Devlin,B.L.,Buenzli,P.R.,Woodruff,M.A.,Simpson,M.J., Newcomputa-\ntionaltoolsandexperimentsrevealhowgeometryaffectstissuegrowthin3Dprintedscaffolds.,\nChemical Engineering Journal, 2023; 475: 145776. doi:10.1016/j.cej.2023.145776.\nLassen, N.E., Andersen, T.L., Pløen G.G., Søe, K., Hauge, E.M., Harving, S., Eschen, G.E.T.,\nDelaisse, J-M.,Coupling of Bone Resorption and Formation in Real Time: New Knowledge\n23\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint \n\nGained From Human Haversian BMUs., Journal of Bone and Mineral Research, 2017; 32:\n1395–1405. doi: 10.1002/jbmr.3091.\nRobling, A., Stout, S.,Morphology of the drifting osteon., Cells Tissues Organs, 1998; 164:\n192–204. doi:10.1159/000016659.\nSchnitzler, C.M., Mesquita, J.M.,Cortical porosity in children is determined by age-dependent\nosteonal morphology., Bone, 2013; 55: 476–486. doi:10.1016/j.bone.2013.03.021.\nCooke, K.M., Mahoney, P., Miszkiewicz, J.J.,Secondary osteon variants and remodeling in\nhuman bone., The Anatomical Record, 2022; 305: 1299–1315. doi:10.1002/ar.24646.\nHegarty-cremer, S.G.D., Borggaard, X.G., Andreasen, C.M., van der Eerden, B.C.J., Simp-\nson, M.J., Andersen, T.L., Buenzli, P.R.,How osteons form: A quantitative hypothesis-\ntesting analysis of cortical pore filling and wall asymmetry., Bone, 2024; 180: 116998. doi:\n10.1016/j.bone.2023.116998.\nTee, S.Y., Bausch, A.R., Janmey, P.A.,The mechanical cell., Current Biology, 2009; 19(17):\nR745 - R748. doi:10.1016/j.cub.2009.06.034.\nMarotti, G., Zallone, A.Z., Ledda, M., Number, Size and Arrangement of Osteoblasts in\nOsteons at Different Stages of Formation., Calcified Tissues, 1976; 21: 96–101. doi:\n10.1007/BF02546434.\nWerner, M., Blanquer, S.B.G., Haimi, S.P., Korus, G., Dunlop, J.W.C., Duda, G.N., Grijpma,\nD.W., Petersen, A., Surface Curvature Differentially Regulates Stem Cell Migration and\nDifferentiation via Altered Attachment Morphology and Nuclear Deformation., Advanced\nscience, 2017; 1600347. doi:10.1002/advs.201600347.\nGudipaty, S.A., Lindbolm, J., loftus, P.D., Redd, M.J., Edes, K., Davey, C.F., Krishnegowda,\nV., Rosenblatt, J.,Mechanical stretch triggers rapid epithelial cell division through piezol.,\nNature, 2017; 118-121. doi:10.1038/nature21407.\nRavasio, A., Cheddadi, I., Chen, T., Pereira, T., Ong, H.T., Bertocchi, C., Brugues, A., Jacinto,\nA., Kabla, A.J., Toyama, Y., Trepat, X., Gov, N., Neves de Almeida, L., Ladoux, B.,Gap\ngeometry dictates epithelial closure efficiency., Nature communications, 2015; 6, 7683. doi:\n10.1038/ncomms8683.\nKuba, S., Simpson, M. J., and Buenzli, P. R. (2026). Stochastic model of curvature-\ncontrolled tissue growth incorporating mechanical cell interactions. GitHub repository.\nhttps://github.com/Shahak-Kuba/2026_Stochastic_Tissue_Growth\n24\n.CC-BY 4.0 International licenseperpetuity. It is made available under a \npreprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in \nThe copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}