Mechanistic interpretation of biological tissue growth experiments with a computational model

–The growth rates of biological tissues are influenced by the existing substrate geometry, mechanobiological processes and the interplay between them. Disentangling the contributions of geometry and mechanobiology experimentally is challenging, as mechanical properties are difficult to measure and tissue samples provide only static snapshot in time. However, the composition of a tissue preserves cues of the dynamic processes that shaped its architecture. In this work, we present a computational model of tissue growth that captures aspects of the interplay between geometry, mechanics, and stochastic biological processes, which we use to generate synthetic tissue compositions directly comparable with experimental samples. This framework enables quantitative analysis of tissue morphology, inference of underlying growth mechanisms, estimation of dynamic rates from single-time-point data, and investigation of how stochasticity contributes to emergent growth patterns. We demonstrate the applicability of the model to simulate the growth of different tissue types by applying this framework to two distinct tissue growth scenarios: (i) tissue grown within 3D-printed porous scaffolds, and (ii) bone formation in cortical pores. Key words:Tissue Engineering, cell culture, bone, single-cell, stochastic 1 Introduction Thegrowthofbiologicaltissuesiscentraltothedevelopment,repair,andregenerationoforgans and tissues. Tissue growth plays an important role for example in bone remodeling, tumor growth,andwoundhealing(Martinetal., 1998;Poujadeetal., 2007;Rollietal., 2012;Wozniak and El Haj,2007; Lowengrub et al.,2009), as well as in tissue engineering and regenerative medicine,wherebiologicaltissuegrowthisinvolvedinengineeringtherepairorreplacementof damaged or diseased tissues (Hollister,2005;O’Brien,2011;Bidan et al.,2013). During tissue ∗Corresponding author. Email address:[email protected] 1 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figure 1– Images from case studies of biological tissue growth. (a)In vitrotissue growth of MC3T3 cells within a400 μm × 400 μm 3D-printed square porous scaffold cultured over 14 days. The extracellular matrix (ECM) is depicted in green and differentiated cell nuclei are indicated in blue. (b)In vivocross section of human bone tissue with multiples cortical pores highlightedinredsurroundingtheHaversiancanalsystemwithembeddedcells(osteocytes)are indicated in black (reproduced with permission fromSkedros et al.(2005)). growth, new biological material is usually created within constrained spaces with experimental observations indicating that the geometry of both the constraining space and the existing tissue substrate influences the rate of growth (Rumpler et al.,2008; Bidan et al.,2012; Schamberger et al.,2023). The accumulation of new material within these geometric constraints induces crowding effects and mechanical stresses on individual cells, which may affect their ability to proliferate, differentiate, and survive through mechanobiological processes (Nelson et al., 2005; Discher et al.,2005; Lim et al.,2010; Ladoux and Mége,2017). Additionally, cell-scale cues arising from the geometry of the existing substrate are also known to influence cell fate (McNamaraetal., 2010;Dobbengaetal., 2016;Callensetal., 2020). Thetissuegrowthdynamics arising from the interplay between geometry, mechanics, and cell fate in turn influence the rate oftissuegrowth,yettheprecisequantitativedetailsoftheseinterplayremainpoorlyunderstood. Experimental approaches for studying the interactions between geometry and mechanics can often be challenging, especially as tissue samples provide single snapshots in time and offer limitedinsightintogrowthdynamics. However,thecompositionofatissuecontainsusefulcues that provide a record of some of the dynamic processes that led to the architecture of a tissue. Biologicaltissuesareprimarilycomposedofanextracellularmatrix(ECM)anddifferentiated cellsembeddedwithinit(Albertsetal., 2002;KurniawanandBouten, 2018). Sincethematerial compositionformsduringgrowth,thepropertiesofdifferentiatedcellsandECMencodevaluable information about the developmental history and can provide a fingerprint of the underlying processes that are regulated in healthy tissues, and dysregulated in diseased tissues. Confocal microscopy and quantitative backscattered electron imaging of tissue samples, such as those in Figure1,provideameanstomeasurematerialproperties,includingspatialinformationrelating to differentiated cells. Several studies use these measurements to characterize ECM properties such as fiber alignment as well as the shape, orientation, and number densities of differentiated cellswithindifferenttissuetypes(Roschgeretal., 2008;Vatsaetal., 2008;Ascenzietal., 2008; van Hove et al.,2009; Britz et al.,2012; Mader et al.,2013; Lanaro et al.,2021). Beyond 2 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint providing a quantitative description of biological tissue composition, these data sets can be leveragedwithinmathematicalmodelingtohelpinferthemechanismsthatgoverntissuegrowth sincethesemodelscanbedesignedtoexplicitlyincorporatespecificmechanisms. Forexample, mineral density distributions in trabecular bone have been used to study bone mineralization dynamics (Ruffoni et al.,2007), osteocyte density in cortical pores to estimate osteoblast burial rates(Buenzli, 2015), andmouseincisorenamelmeasurementstoformulatepredictivetheories ofcellularresponsestostrainvariations(Coxetal., 2024). Inthisregard,mathematicalmodeling offersapowerfulframeworkforgeneratingandtestinghypotheses,analyzingexperimentaldata, and could be useful for disentangling interplay between geometry, mechanics, and cell fate in the processes underlying tissue growth. Several mathematical models of tissue growth have been developed to predict and help understand the evolution of biological tissues, with varying levels of biological detail. Phe- nomenologicalmodelscapturethemacroscopicbehaviorofgrowingtissuesbydefininggrowth lawsthatlinkthevelocityofthetissueboundaryanditscurvature,forexample,meancurvature flow (Rumpler et al.,2008; Guyot et al.,2014; Callens et al.,2023). Continuum mechan- ics approaches typically employ growth tensors to describe tissue deformation and expansion (Ambrosi and Guana,2007; Dunlop et al.,2010; Goriely, 2007). However, these growth laws and growth tensors are difficult to relate to biological processes. Cell-based models, which include both continuum and discrete approaches, complement phenomenological descriptions by explicitly accounting for cell populations and modeling their associated behaviors (Bidan et al.,2012; Joly et al.,2013; Gahffarizadeh et al.,2018). These modeling frameworks have shown that some curvature effects arise from the crowding and spreading of cells (Alias and Buenzli,2017;Hegarty-Cremeretal., 2021;Kubaetal., 2026). Attheindividualcelllevel,these models are also able to connect mechanical stresses to tissue organization and to cellular events (Baker et al.,2019; Murphy et al.,2019; Buenzli et al.,2025; Brown et al.,2025). While these modeling approaches have provided valuable insights into tissue growth, they tend to focus on theevolutionofthetissuesurfaceratherthangeneratingatissuecompositionthatiscomparable to experimental samples. Inthiswork,wepresentacomputationalmodelofbiologicaltissuegrowththatcanbeusedto simulateandanalyzethedynamicprocessesthatleadtoobservedtissuecompositionsinconfined geometries. The model explicitly represents individual cells and includes mechanical interac- tions between cells, the generation of ECM, and incorporates cell proliferation, differentiation, and death as stochastic events. Unlike previous interface-based or curvature-controlled growth models, this framework explicitly records the positions and orientations of differentiated cells as tissue forms, enabling direct comparison of synthetic and experimental tissue compositions at single-cell resolution. Comparisons between model-generated predictions and experimental dataenablesusto: (i)relatemechanicalinteractionsbetweencellstoboththeirorientationsand the shape of the tissue interface; (ii) identify key mechanisms that are important for regulating tissue growth; (iii) estimate cell differentiation rates at different times from single snapshot; 3 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint and,(iv)analyzehowmodeledmechanismsandstochasticityshapetissuecompositionandgive rise to random growth patterns. To demonstrate the flexibility of this framework, we apply it toinvestigatedatafromtissueengineeringcellcultureexperimentsin3D-printedscaffolds,and bone formation in cortical pores. 2 Materials & Methods This study is primarily simulation-based, using comparisons between computational model predictionsandexperimentalobservationstointerpretthemechanismsunderlyingtissuegrowth. While experimental data provide quantitative description of tissue composition, they do not reveal the precise mechanisms that control the growth behaviour. By applying mechanistic simulation models to interpret these data, we obtain biological insights that are otherwise difficult to access. 2.1 Mathematical model The mathematical model simulates tissue growth by evolving a discrete representation of the tissueinterface,withnewlyformedtissuefillingtheregionbehindit. Theinterfaceisrepresented as a one-dimensional chain of𝑁 tissue-forming cells connected through their boundaries. The positionofthe 𝑖thcellattime 𝑡isgivenbythetwo-dimensionalpositionvectorsofitsboundaries, 𝒓𝑖−1(𝑡) and 𝒓𝑖 (𝑡) (Figure2a). Tissuegrowtharisesfromthreeprocesses: (i)mechanicalinteractionsbetweenneighbouring cells, (ii) ECM deposition at the interface, and (iii) stochastic cellular events. Following previous discrete cell–cell mechanical interaction modeling approaches, intercellular forces are represented by representing cell bodies as mechanical springs (Murray et al.,2009; Murphy et al.,2019; Buenzli et al.,2025). Cells are assumed to generate an area of tissue material at a constant rate𝑘f (μm2/day) that displaces them perpendicularly to the interface (Buenzli, 2015; Alias and Buenzli,2017; Kuba et al.,2026). Together, these two mechanisms yield two independentvelocitycontributionsthatevolvethepositionofagivencellboundary 𝒓𝑖,according to d𝒓𝑖 d𝑡 =𝒗 m 𝑖 +𝒗 f 𝑖, 𝑖=1,2, ..., 𝑁,(1) where𝒗 m 𝑖 represents tangential motion arising from mechanical relaxation, and𝒗f 𝑖 represents normal motion due to tissue formation. Thevelocity 𝒗m 𝑖 isconstructedundertheassumptionofoverdampedmotion,appropriatefor slow cellular dynamics in viscous environments (Purcell,1977). In this regime, inertial effects areneglectedandvelocityisproportionaltothenetforceactingontheboundary 𝒓𝑖 (𝑡) suchthat 𝒗m 𝑖 = 1 𝜂 𝑭(net) 𝑖 ,(2) 4 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figure2 –Overviewofthemathematicalmodel.(a)Tissuegrowthregion R andtime-dependent tissue area 𝛺(𝑡) showing differentiated cells (red) and the discrete cell-based tissue interface 𝑆(𝑡). Tissue formation displaces interface cells perpendicularly in the direction ¯𝒏𝑖 with a velocityproportionaltothelocalcelldensity 𝑞𝑖. (b)Mechanicalforcediagramatcellboundary 𝒓𝑖 showing two restoring forces from springs on either side and the normal reaction force as well as the net force. (c-f) Examples of discrete cells undergoing proliferation, asymmetric proliferation, apoptosis, and differentiation over a time stepΔ𝑡. Highlighted cells undergo the indicated process, and red cells denote differentiated cells as shown in (a). where 𝜂 is the viscous drag coefficient. The net force is composed of three contributions (Figure 2b): two restoring forces,𝑭(+) 𝑖 and 𝑭(−) 𝑖 , from the cells on either side, and a reaction force 𝑭(n) 𝑖 from the existing substrate. The restoring forces act along the cell bodies, with magnitudesdeterminedbyaprescribedforcelaw 𝑓 (ℓ𝑖 (𝑡)) whichmaybeHookeanornonlinear, and ℓ𝑖 (𝑡) = ∥ 𝒓𝑖 (𝑡) − 𝒓𝑖−1(𝑡) ∥ is the cell length (Murray et al.,2012; Murphy et al.,2019). he reaction force is defined to oppose the normal component of the restoring forces, thereby restricting mechanical relaxation to the tangent direction of the interfaceˆ𝝉𝑖. Consequently, the net force is given by the tangential projection 𝑭(net) 𝑖 = h 𝑭(+) 𝑖 + 𝑭(−) 𝑖  · ˆ𝝉𝑖 i ˆ𝝉𝑖. (3) Thevelocitycontribution 𝒗f 𝑖,modelstheperpendiculardisplacementofcellsinthedirection ¯𝒏𝑖 as they deposit ECM at the interface (Figure2a). Compatibility considerations requires that 5 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint the magnitude of this velocity contribution to satisfy ∥𝒗f 𝑖 ∥=𝑘 f 𝑞𝑖,(4) where𝑞 𝑖 =1/ℓ 𝑖 is the density of the𝑖th cell (Buenzli,2015; Alias and Buenzli,2017; Hegarty- Cremeretal., 2021;Kubaetal., 2026). Inthisdiscreteframework,perpendiculardisplacements of the cell segments can cause them to overlap or separate, which leads to a disconnected interface. To combat this issue we implement a predictor-corrector algorithm to calculate𝒗f 𝑖 depending on whether neighboring segments intersect or separate after being displaced. A full description of this algorithm is described in detail inKuba et al.(2026). Inadditiontothemechanicalandformationprocesses,themodelincorporatesfourstochastic cell events: • Proliferation(Symmetricdivision): acelldividesintotwoactivetissue-formingdaugh- ter cells that both lie on the tissue interface (Figure2c). • Proliferation (Asymmetric division): a cell divides into two daughter cells where one remains an active tissue-forming cell and the other differentiates and embeds into the ECM (Figure2d). • Apoptosis: a cell dies and is then removed from the tissue interface (Figure2e). • Differentiation: a cell differentiates and embeds into the ECM (Figure2f). A key novelty of the framework is that it explicitly records the positions of differentiated cells (marked in red in Figure2a,d,f) during asymmetric division and differentiation, enabling single-cell resolution comparisons between model predictions and experimental observations. Although asymmetric division and direct differentiation both generate differentiated cells, they represent distinct cellular mechanisms. While these processes are not easily distinguishable in static experimental observations, they can be independently controlled and analysed within the computational framework. These events are assumed to occur over a sufficiently small time intervalΔ𝑡, such that at most one event takes place within[𝑡, 𝑡 + Δ𝑡) (Baker et al.,2019; Murphy et al.,2020). Each cell may undergo proliferation (either symmetric division or asymmetric division), apoptosis, or differentiation with constant rates (probabilities per unit time)𝑃sym, 𝑃asym, 𝐴, and 𝐷, re- spectively. Stochastic implementation proceeds by generating three random numbers from a uniformdistributiontodetermine: (i)whetheraneventoccurs, (ii)whicheventoccurs, and(iii) which cell undergoes the event (Gelman et al.,2013; Baker et al.,2019; Murphy et al.,2020). Implementation details are provided inA.1 and open-source code is available atKuba et al. (2026). 6 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint 2.2 Connecting model parameters to experimental observations To enable comparison between model-generated and experimental tissue compositions, model parameters were calibrated to reproduce experimentally measured differentiated cell densities. In the model, differentiated cells arise at the moving tissue interface through either direct differentiation (𝐷) or asymmetric division (𝑃asym), such that the total differentiation rate is 𝐷tot = 𝐷 + 𝑃asym. (5) FollowingBuenzli(2015), experimentally consistent differentiated cell densities𝜌 are obtained by approximating 𝐷tot = 𝑘f 𝜌, (6) where 𝑘f is the tissue formation rate. For each case study, experimentally estimated values of𝜌 and 𝑘f (Section 2.3) were substituted into Eq. (6) to determine𝐷tot. To characterize the relative contribution of differentiation pathways, we define 𝛼 = 𝐷 𝐷tot ,(7) representing the fraction of differentiation occurring through direct differentiation. Since direct differentiationremovesanactivetissue-formingcellfromtheinterface,increasing𝛼reducesthe density of active cells and influences interface morphology via the normal velocity formulation (Eq. (4)). The evolution of the active tissue-forming cell population𝑁is described by d𝑁 d𝑡 =𝛽𝑁,(8) where𝛽=𝑃 sym −𝐴−𝐷 is the net growth rate. Here,𝛽 > 0 corresponds to a proliferative population and𝛽 < 0 to a declining population. Together,𝛼 and 𝛽 provide independent control over differentiation mode and population dynamics in the simulations 2.3 Experimental data Previouslypublishedexperimentalmeasurementsofdifferentiatedcelldensityandextracellular matrix (ECM) formation rates were used to calibrate model parameters and ensure consistency between simulated and observed tissue compositions. Two experimental contexts were con- sidered: (i) tissue engineering cultures within 3D-printed constructs, and (ii) bone formation within cortical pores. For tissue engineering applications, data were drawn from in vitro studies of MC3T3-E1 osteoblasticcellsculturedinthin3D-printedporousscaffolds(Buenzlietal., 2020;Lanaroetal., 2021; Browning et al.,2021). In those experiments, cells were seeded within square pores of 7 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint side length200–600μm and cultured for up to 28 days. Tissue growth was monitored using DAPI staining of cell nuclei and actin filaments, with ethidium homodimer staining indicating theabsenceofcelldeath(Buenzlietal., 2020). Imageswereacquiredatmultipletimepoints(4, 7, 10, 14, 18, and 28 days), from which differentiated cell density, cell orientation, tissue area coverage, and pore infill duration were quantified (Lanaro et al.,2021; Browning et al.,2021). Althoughorientation,coverage,andinfilltimevariedwithporesize,parameterestimationacross allporesindicatedthatatlatetime( 𝑡 > 14days)thedensityofdifferentiatedcellsremainaround 𝜌 ≈ 0.003 cells/μm2 (Browning et al.,2021). The corresponding ECM formation rate was estimated as𝑘f ≈ 87.84 μm2/day (Kuba et al.,2026). For cortical bone formation, matrix apposition rate (MAR) and osteocyte density (𝜌) were obtainedfrompreviouslypublisheddouble-labelingexperimentsandimagingstudiesinanimal models and human bone (Metz et al.,2003; Hannah et al.,2010; Power et al.,2012). Reported MAR values in human cortical bone range between0.5–3μm/day (Jones,1974; Martin et al., 1998; Buenzli et al.,2014; Alias and Buenzli,2018). To ensure consistency with this range (target ≈ 1.5 μm/day), the model parameter𝑘f was set to30 μm2/day, assuming a typical osteoblast size of∼ 20 μm and using Eq. (4). Consistent with reported estimates (Alias and Buenzli, 2018; Parfitt, 1983; Hannah et al.,2010; Power et al.,2012), osteocyte density was taken as𝜌 = 0.000625 cells/μm2 throughout the pore. 3 Results We apply the computational model to two case studies: (i) tissue growth within 3D-printed scaffolds (Section3.1) and (ii) bone formation within cortical pores (Section3.2). In the tissue engineering study, we investigate how differentiation mechanisms, cell population dynamics, and mechanical properties of cells shape tissue growth patterns, interface morphology, and differentiated cell orientations. In the bone formation case study we examine how stochasticity influences the distribution of differentiated cells and the structural symmetry of the pore. 3.1 Tissue Engineering To reflect the experimental setup described in Section2.3, the tissue interface was initialized in a 300 μm × 300 μm square with 𝑁 = 60 uniformly distributed cells of length20 μm, approximatingosteoblasticsize. FromEq.(6),weestimated 𝐷tot ≈ 0.26,andduetotheabsence of cell death we set𝐴 = 0. Since the fraction of differentiated cells (𝛼) that emerge from directdifferentiationandthenetgrowthrateofthetissue-formingpopulation( 𝛽)aredifficultto measure experimentally, we investigated their effects by varying these parameters. Figure 3a shows representative model-generated tissue compositions, with differentiated cells shown in red and the tissue interface indicated at the initial time (dark blue) and at 90% pore fill (blue). Within the newly formed tissue, differentiated cells vary in orientation and size 8 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figure 3– Example realisations of the tissue growth model for a series of control parameters 𝛼and𝛽values. Simulations are initialised with a square geometry to emulate the cell culture experiments with𝐷 tot =0. 003. (a) Realisations of the model for different values of𝛼 with 𝛽 = 0. (b)Exampleimagesofexperimentsculturedwithing 300 μm × 300 μm,reproducedwith permissions fromBuenzli et al.(2020). (c) Solid lines are the average normalized area of void space over time for different values of𝛽 with 𝛼 = 0.2, and black points are averages estimated fromexperimentsculturedinscaffoldsofthesamesize. Thesimulatedaveragesweregenerated from from a sample size of 50 realizations for each value of𝛽. 9 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figure4–Simulationsofthestochasticmodelshowingtheeffectofcellstiffnessonembedded cell orientations. Each column shows a single realization for (a) low (𝑘=5 N/μm) (b) intermediate (𝑘 = 100 N/ μm) and (c) high (𝑘 = 500 N/ μm) stiffness. The top row displays tissue growth with differentiated cells (red) and the interface at𝑡 = 0, 5, and 10 days, using a color gradient from dark blue (early) to light blue (late). Subsequent rows show the frequency of cell orientations measured relative to the x-axis and binned into 24 intervals over[0◦ , 360◦ ). 10 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figure 5– Comparison of experimental and simulated cell orientation distributions. (a) Mea- suredcellorientationsfromcellcultureexperimentsondays 4, 7, and 10ina 300 μm × 300 μm 3D-printed square porous scaffold reported byLanaro et al. (2021) (reproduced with per- missions). (b-d) Cumulative cell orientation frequencies from 100 realizations for low (𝑘 = 5 N/ μm), intermediate (𝑘 = 100 N/ μm) and high (𝑘 = 500 N/ μm) cell stiffness at simulation times that correspond to the experimental data. butremainapproximatelyuniformlydistributed,despitetheincreasingdensityoftissue-forming cells as the interface evolves. To investigate the influence of differentiation modes, in Figure3a wevaried 𝛼andset 𝛽 = 0sothat 𝑃sym, 𝑃asym,and 𝐷 balancetomaintainaconsistentpopulation of tissue-forming cells. Comparison with experimental observations (Figure3b) indicates that 𝛼 = 0.2 provides the closest agreement between model predictions and experiments, whereas larger 𝛼 values yield a rougher and less consistent interface morphology. Figure3c shows the effectofadynamictissue-formingcellpopulationonthetimerequiredtoachieveporeinfilling. In these simulations,𝛼 was fixed at0.2, and growing (𝛽 > 0), constant (𝛽 = 0) and shrinking (𝛽 < 0) populations were considered. As in Figure3a, simulations were performed until 90% pore fill, and the normalized tissue area, averaged from 50 realizations, was plotted over time for each𝛽 value alongside the averaged experimental data for scaffolds of this size. Figure4examinestheinfluenceofcellmechanicalstiffnessondifferentiatedcellorientation. Stiffnesslevelswerecategorizedaslow( 𝑘 = 5 N/μm),intermediate( 𝑘 = 100 N/μm),andhigh (𝑘 = 500 N/μm). Representativerealizationsandcorrespondingroseplotsshowthatorientation distributions are initially anisotropic across all stiffness levels (𝑡 = 4 days), but anisotropy de- 11 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint creasesovertimeforlowerstiffnessvalues. Comparisonofcumulativeorientationdistributions with experimental measurements (Lanaro et al., 2021) indicates improved agreement at lower stiffness, although the associated interface morphologies do not reproduce the experimentally observed rounding (Figure 5). 3.2 Bone formation To model bone tissue formation within irregular resorption cavities, the tissue interface was initialized as a deterministically perturbed circle. Specifically,𝑁values were sampled from a uniform distribution and smoothed with respect to the angular coordinate𝜃 ∈ [ 0, 2𝜋) using the loessfunction to generate spatially correlated perturbations (Cleveland and Grosse,1991). The resulting smoothed values𝑢𝑖 ∈ [ 0, 1] were scaled by a perturbation amplitude˜𝑅 and added to a circle of radius𝑅0 such that, at𝑡 = 0, 𝑅𝑖 = 𝑅0 + 𝑢𝑖 ˜𝑅, 𝑖 = 1, . . . , 𝑁 . (9) Experimental evidence indicates that bone formation initiates once osteoblast density reaches a threshold, at which point proliferation ceases (Lassen et al.,2017). Accordingly, proliferation rates were set to zero (𝑃sym = 𝑃asym = 0), so differentiation occurred solely through direct differentiation, yielding𝐷 = 𝑘f𝜌 ≈ 0.019. The apoptosis rate𝐴 was chosen to be lower than𝐷 to prevent rapid cell loss. Figure 6b shows a representative realization for an artificial resorption cavity, with osteo- cytes marked in black and growth regions shaded according to formation time, consistent with experimental images (Figure6a). Osteocyte density was quantified within regions formed over a time interval of𝑇 = 7 days, corresponding to lamellar deposition in cortical bone. Summary statistics (Figure6c) demonstrate that local osteocyte density remains consistent and closely matches the prescribed density𝜌, despite variation in region shape and size over time. AlthoughtherealizationinFigure 6bshowssymmetricporeclosure,corticalporesmayalso developasymmetricstructures. Asymmetrywasquantifiedusingthewallthicknessasymmetry ratio W.Thasym, defined using the pair of opposing wall thicknesses with the greatest disparity (Figure6a), where wall thickness is measured as the distance between the inner and outer pore boundaries (Hegarty-Cremer et al.,2024): W.Thasym = W.Thmin W.Thmax .(10) Figure 6d shows representative realizations classified by W.Thasym, and Figure 6e summarizes osteocyte densities in the wall regions used to compute this ratio. These results indicate that osteocyte density is higher near thinner walls, and the most frequently observed asymmetries fall within the range W.Thasym ∈ [0.6,0.7). 12 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figur e 6– Model simulations of bone formation in an artificial osteon. (a) Experimental image illustrating the measurement of W.Thmax and W.Thmin used to compute the W.Thasym ratio. (b) Representative model realization in an artificial osteon, with osteocytes marked in black. Pseudo-lamellae formed over a time interval𝑇 = 7 days are shaded according to formation time. (c) Box plots of embedded cell densities within wall regions formed over 𝑇 days from 500 realizations. Shading corresponds to the wall regions illustrated in (a). (d) Example model realizations spanning a range of W.Thasym values. (e) Summary statistics of embedded cell densities measured in rectangular regions nearW.Thmax (blue) andW.Thmin (red), grouped by 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 and total differentiation rate is𝐷tot = 𝑘f 𝜌 where 𝜌 = 0.000625 cells/μm2. 13 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint 4 Discussion Experimental evidence indicates that biological tissue growth is strongly regulated by substrate geometry and mechanobiological cues (Nelson et al.,2005; Rumpler et al.,2008; Bidan et al., 2012;LadouxandMége, 2017;Callensetal., 2020). Inconfinedgeometries,spatialconstraints causetissue-formingcellstocrowdorspreadduringgrowth,generatingmechanicalstressesthat inturninfluencecellularprocessessuchasproliferationanddifferentiation(Nelsonetal., 2005; Werner et al.,2017; Gudipaty et al.,2017). However, disentangling the mechanistic interplay betweengeometry,mechanics,andcellularprocessesintissuegrowthremainschallengingwhen relyingsolelyonexperimentaltissuesamples,eventhoughsuchsamplescontainfingerprintsof theprocessesthatshapedthem. Toaddressthislimitation,weproposeacomputationalmodelof tissuegrowththatgeneratessynthetictissueswithsingle-cellresolution,enablingdirectcontrol of governing mechanisms and systematic comparison with experimental data to isolate their individual contributions. Duringtissuegrowth,cellsmayundergoproliferation,differentiation,andapoptosis,withthe relative contributions of these processes varying across tissue types. Histological staining can reveal the presence or absence of certain processes, but it cannot resolve the specific processes that occurred during growth. In tissue engineering experiments (Figure1a), DAPI staining indicates that there is a substantial amount of cell proliferation during growth, while ethidium homodimer staining indicates the absence of cell death. Model simulations indicate that, given the incorporated mechanisms, closer agreement with experimental observations and smoother, more consistent tissue interfaces are achieved when the majority of differentiated cells arise from asymmetric division (𝛼 = 0.2) rather than direct differentiation (Figure3a). In addition, comparisonofnormalizedtissueareaovertimeindicatesthatanincreasingpopulationoftissue- forming cells more accurately captures late-stage growth dynamics (Figure3c), consistent with previous quantitative and modeling studies (Buenzli et al.,2020; Browning et al.,2021). Inthesetissueengineeredexperiments,differentiatedcellsareassumedtobecomequiescent onceembeddedwithintheECM.Consequently,differentiatedcellorientationsprovidearecord of the position and shape of the tissue interface at the time of differentiation. Previous mod- eling studies of evolving tissue interfaces have demonstrated a relationship between interface morphologyandthestrengthofmechanicalinteractionsbetweencells(AliasandBuenzli, 2017; Kubaetal., 2026). Consistentwiththis, modelsimulationsshowthatintermediatecellstiffness produces tissue interfaces that more closely resemble the experimentally observed rounding (Figure 4), whereas low stiffness yields orientation distributions with reduced anisotropy that better match experimental measurements (Lanaro et al.,2021) (Figure5). This discrepancy suggests that additional mechanisms may contribute to tissue growth mechanics beyond those captured by cell–cell mechanical interactions alone. A plausible contributor to these additional mechanisms is surface tension. In their experi- mental study,Lanaro et al.(2021) observed groups of cells spanning perpendicular fibers and 14 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Figure 7– Example realisations of the tissue growth model in different initial geometries. The geometries are motivated by previously conducted cell culture experiments ofRumpler et al. (2008) (triangle),Bidan et al.(2013) (cross), andVandenHeuvel et al.(2023) (wavy). proposed that cells actively seek nearby attachment points to maintain surface tension during growth. This behavior is consistent with purse–string–based closure mechanisms observed in wound healing (Ravasio et al.,2015) and with evidence that surface tension plays a central roleintissueorganizationandmechanicalregulation(Bidanetal., 2016). Inthecurrentmodel, surfacetension–relatedforcesactingperpendiculartotheinterfacearenotrepresentedduetothe definitionofthereactionforce 𝑭(n) 𝑖 . Thesefindingsunderscoretheimportanceofincorporating surface tension–mediated effects into tissue growth models, as neglecting them may obscure keymechanicalcontributionstointerfaceevolution. Althoughthiscasestudyfocusedonsquare geometries,thesameanalysiscanbeextendedtootherporeshapes(Figure 7)tofurtherexamine the role of surface tension in tissue growth. InSection 3.2, weinvestigatedhowstochasticityinfluencesthespatialdistributionofdiffer- entiated cells and the resulting growth patterns during bone formation within cortical pores. In this setting, a finite population of osteoblasts deposits new matrix at the pore surface without proliferation, differentiating into osteocytes that become embedded within the mineralized tis- sue. Experimental studies report that osteocyte density decreases with distance from the initial pore boundary (Power et al.,2012), although differentiation rates cannot be directly measured. Consistent with previous modeling work (Buenzli,2015), our simulations show that osteocyte density within newly formed lamellar regions remains consistent with the prescribed constant 15 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint value𝜌, despite variations in regional shape and size (Figure 6b,c). This suggests that dy- namic differentiation rates could be inferred from single experimental snapshots by quantifying osteocyte densities within lamellae and rearranging Eq. (6), providing a basis for incorporat- ing curvature-, stress-, or property-dependent differentiation rules (e.g.,𝐷 (𝜅, 𝜎, . . .)) into the modeling framework. CorticalboneformationoccursaroundHaversiancanalswhosepositionsareoftenasymmet- ric(Robling, 1998;SchnitzlerandMesquita ,2013;Cookeetal., 2022),motivatinginvestigation into the mechanisms underlying such asymmetries. Previous modeling studies have attributed asymmetricgrowthtobiologicalheterogeneities, suchasspatialvariationsinmatrixdeposition rates (Hegarty-Cremer et al.,2024). In contrast, our analysis (Figure6d,e) demonstrates that asymmetriescanarisepurelyfromstochasticdifferentiationandapoptosisofbone-formingcells. Thediscreteformulationofthemodelcapturesdensityfluctuationsassociatedwiththeseevents, leading to locally reduced normal velocities (Eq. (4)) and emergent interface irregularities. These findings suggest that at least part of the variability observed in cortical pore morphology mayresultfromintrinsic stochasticcellularprocesses, althoughfurtherexperimentalvalidation is required. 5 Conclusions and future work Experimental evidence demonstrates that tissue growth is governed by the interplay between geometry,mechanics,andbiologicalprocesses(Nelsonetal., 2005;Rumpleretal., 2008;Bidan et al.,2012; Ladoux and Mége,2017; Callens et al.,2020). Disentangling these mechanisms experimentally remains challenging, as mechanical properties are difficult to measure (Tee et al.,2009) and tissue samples provide only static snapshots of a dynamic process. However, the material deposited behind the advancing tissue interface retains structural signatures of the processes that shaped tissue architecture. Here, we introduced a computational model thatsimulatestissuegrowthatsingle-cellresolutionandgeneratessynthetictissuecompositions directlycomparablewithexperimentalobservations. Throughtwocasestudies,wedemonstrated thattheframeworkcan(i)identifythecontributionofspecificmechanismstotissuemorphology, (ii)estimatedynamicratesfromsingle-time-pointdata,and(iii)generatetestablehypothesesfor the origins of stochastic growth patterns. These results highlight the utility of model-generated tissue compositions as a complementary tool for investigating tissue growth dynamics and enabling in silico hypothesis testing. Futureextensionscouldfocusontailoringtheframeworktospecifictissuetypesbyincorpo- ratingbiologicallyinformedprocessrules. Forexample,inboneformation,differentiationrates could be extended to depend on interface curvature or mechanical stress and calibrated using experimental estimates. Similar extensions could be applied to proliferation, apoptosis, and the tissue formation rate𝑘f, consistent with previous studies in which𝑘f depended on interface radiusorporosity(AliasandBuenzli, 2018), andcelleventratesdependedoncelllength(Mur- 16 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint phy et al.,2020). Parameter estimation could then be performed by comparing synthetic and experimental datasets using statistical inference techniques. Additionally, the model’s ability to track embedded cell positions enables spatial statistical analyses, such as point process or pair-correlationmethods,providingaquantitativeframeworkfortestingcompetingmechanistic hypotheses and guiding future experimental design.

ThisresearchwassupportedbytheAustralianResearchCouncil(DP190102545,DP230100025). SK and PRB acknowledge support from the Max Planck Queensland Centre for the Materials ScienceofExtracellularMatrices. WethankBrennaDevlinandMariaA.Woodruffforproviding the experimental image presented in Figure 1a. Appendix A Mathematical Model A.1 Implementation of stochastic cellular processes in mathematical model To simplify notation, we write the probability that the𝑖th cell proliferates is𝑃sym 𝑖 (symmetric division). Similarly for asymmetric proliferation, apoptosis, and differentiation, we write the probabilities as 𝑃asym 𝑖, 𝐴𝑖, and 𝐷𝑖. The algorithm at each numerical time stepΔ𝑡 has been implemented as follows: 1. Generate 𝑝1, 𝑝2, 𝑝3 from a uniform distribution in[0, 1]. 2. If 𝑝1 ≤ Í𝑁 𝑖=1(𝑃sym 𝑖 + 𝑃asym 𝑖 + 𝐴𝑖 + 𝐷𝑖)Δ𝑡, where𝑁 isthetotalnumberofcellsthatform the interface, then one of the cells is experiencing one of the cellular processes over the time interval[𝑡 + Δ𝑡). Otherwise continue to the next simulation time step. 3. Using 𝑝2, we determine which cellular process occurred under the following conditions: • If 𝑝2 ≤ Í𝑁 𝑖=1 𝑃sym 𝑖 Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) then the occurring cellular process is proliferation. •Else if Í𝑁 𝑖=1 𝑃sym𝑖 Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) < 𝑝 2 ≤ Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 )Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) then the occurring cellular process is asymmetric proliferation. •Else if Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 )Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) < 𝑝 2 ≤ Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖)Í𝑁 𝑖=1 (𝑃sym𝑖 +𝑃asym𝑖 +𝐴 𝑖+𝐷 𝑖) then the occurring cellular process is apoptosis. •Else the occurring cellular process is differentiation. 4. Using𝑝 3,wedeterminewhichcellexperiencedthecellularprocess. Forexample,assum- ingthatthepreviousstepresultedinproliferation(symmetricdivision),weassertthatthe 17 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint 𝑗th cell proliferates when the conditions below are met, 0< 𝑝 3 ≤ Í 𝑗 𝑖=1 𝑃sym𝑖 Í𝑁 𝑖=1 𝑃sym𝑖 , 𝑗=1, Í 𝑗−1 𝑖=1 𝑃sym𝑖 Í𝑁 𝑖=1 𝑃sym𝑖 < 𝑝 3 ≤ Í 𝑗 𝑖=1 𝑃sym𝑖 Í𝑁 𝑖=1 𝑃sym𝑖 ,2≤ 𝑗 ≤ 𝑁 . A similar process is followed in the case that the previous step results in apoptosis, embedment, and simultaneous proliferation and embedment. 5. Based on cell event, we either add or remove cell boundaries from the interface, linearly interpolating the boundary nodes over the region where the𝑗th cell spanned as described in Section2.1.

Poujade, M., Grasland-Mongrain, E., Hertzog, A., Silberzan, P.,Collective migration of an epithelialmonolayerinresponsetoamodelwound.,CellBiology,2007;104: 15988–15993. doi: 10.1073/pnas.0705062104. Rolli,C.G.,Nakayama,H.,Yamaguchi,K.,Spatz,J.P.,Kemkemer,R.,Nakanishi,J., Switchable adhesive substrates: Revealing geometry dependence in collective cell behavior., Biomateri- als, 2012; 33: 2409–2418. doi:10.1016/j.biomaterials.2011.12.012. Martin, R.B., Burr, D.B., Sharkey, N.A., Fyhrie, D.P.Skeletal Tissue Mechanics., New York: springer, 1998. Lowengrub,J.S.,Frieboes,H.B.,Jin,F.,Chuang,Y-L.,Li,X.,Macklin,P.,Wise,S.M.,Cristini, V.,Nonlinearmodellingofcancer: bridgingthegapbetweencellsandtumours.,Nonlinearity, 2009; 23: R1. doi:10.1088/0951-7715/23/1/R01 . Wozniak,P.,Elhaj,A.J., TissueEngineeringUsingCeramicsandPolymers,Chapter14-Bone regeneration and repair using tissue engineering., Woodhead Publishing, 2007. Hollister, S.J.,Porous scaffold design for tissue engineering, Nature Materials, 2005; 4: 518– 524. doi:10.1038/nmat1421. O’Brien, F.J.,Biomaterials & scaffolds for tissue engineering, Materials Today, 2011; 14: 88–95. doi:10.1016/S1369-7021(11)70058-X. Bidan,C.M.,Kommareddy,K.P.,Rumpler,M.,Kollmannsberger,P.,Fratzl,P.,Dunlop,J.W.C., GeometryasaFactorforTissueGrowth.,AdvancedHealthcareMaterials,2013;2: 186–194. doi: 10.1002/adhm.201200159. 18 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Rumpler, M., Woesz, A., Dunlop, J.W.C, van Dongen, J.T, Fratzl, P.,The effect of geometry on three-dimensional tissue growth., Journal of The Royal Society Interface, 2008; 27: 1173– 1180. doi: 10.1098/rsif.2008.0064. Bidan, C., Kommareddy, K., Rumpler, M., Kollmannsberger, P., Bréchet, Y., Fratzl, P., Dunlop, J.,How Linear Tension Converts to Curvature: Geometric Control of Bone Tissue Growth., PloS one, 2012; 7: e36336. doi: 10.1371/journal.pone.0036336. Schamberger, B., Ziege, R., Anselme, K., Ben Amar, M., Bykowski, M., Castro, A.P.G., Cipitria, A., Coles, R.A., Dimova, R., Eder, M., Ehrig, S., Escudero, L.M., Evans, M.E., Fernandes, P.R., Fratzl, P., Geris, L., Gierlinger, N., Hannezo, E., Iglič, A., Kirkensgaard, J.J.K.,Kollmannsberger,P.,Kowalewska,Ł.,Kurniawan,N.A.,Papantoniou,I.,Pieuchot,L., Pires, E.F., Roschger, A., Bidan, C.M., Dunlop, J.W.C.,Curvature in Biological Systems: Its Quantification, Emergence, and Implications across the Scales., Advanced Materials, 2023; 35: 2206110. doi:10.1002/adma.202206110. Callens,S.J.P.,Uyttendaele,R.J.C.,Fratila-Apachitei,L.E.,Zadpoor,A.A., Substratecurvature as a cue to guide spatiotemporal cell and tissue organization., Biomaterials, 2020; 232: 119739. doi:10.1016/j.biomaterials.2019.119739. McNamara,L.E.,McMurray,R.J.,Dalby,M.J., NanotopographicalControlofStemCellDiffer- entiation., Journal of Tissue Engineering, 2010; 1: 120623. doi:10.4061/2010/120623. Dobbenga, S., Fratila-Apachitei, L.E., Zadpoor, A.A.,Nanopattern-induced osteogenic differ- entiation of stem cells – A systematic review., Acta Biomaterialia, 2016; 46: 3–14. doi: 10.1016/j.actbio.2016.09.031. Nelson, C.M., Jean, R.P., Tan, J.L., Liu, W.F., Sniadecki, N.J., Spector, A.A., Chen, C.S., Emergentpatternsofgrowthcontrolledbymulticellularformandmechanics.,Proceedingsof theNationalAcademyofSciences,2005;102: 11594–11599.doi: 10.1073/pnas.0502575102. Discher, D.E., Janmey, P., Wang, Y-L.,Tissue Cells Feel and Respond to the Stiffness of Their Substrate., American Association for the Advancement of Science, 2005; 310: 1139–1143. doi: 10.1126/science.1116995. Ladoux,B.,Mége,R.M., Mechanobiologyofcollectivecellbehaviours.,NatureReviewsMolec- ular Cell Biology, 2017; 18: 743–757. doi:10.1038/nrm.2017.98. Lim, C.T., Bershadsky, A., Sheetz M.P.,Mechanobiology., Journal of The Royal Society Inter- face, 2010; 7: S291–S293. doi:10.1098/rsif.2010.0150.focus. Alberts, B., Johnson, A., Lewis, J., Raff, M., Roberts, K., Walter, P.,Molecular Biology of the Cell, 4th edition., New York: Garland Science, 2002. 19 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Kurniawan, N.A., Bouten, C.V.C.,Mechanobiology of the cell–matrix interplay: Catching a glimpseofcomplexityviaminimalisticmodels,ExtremeMechanicsLetters,2018;20: 59–64. doi: 10.1016/j.eml.2018.01.004. Skedros,J.G.,Holmes,J.L.,Vajda,E.G.,Bloebaum,R.D., Cementlinesofsecondaryosteonsin humanbonearenotmineral-deficient: Newdatainahistoricalperspective.,TheAnatomical Record, 2005; 286: 781–803. doi:10.1002/ar.a.20214. Roschger,P.,Paschalis,E.P.,Fratzl,P.,Klaushofer,K., Bonemineralizationdensitydistribution in health and disease., Bone, 2008; 42: 456–466. doi:10.1016/j.bone.2007.10.021. Vatsa, A., Breuls, R.G., Semiens, C.M., Salmon, P.L., Smit, T.H., Klein-Nulend, J.,Osteocyte morphology in fibula and calvaria — Is there a role for mechanosensing?, Bone, 2008; 43: 452–458. doi:10.1016/j.bone.2008.01.030. Ascenzi, M., Gills, J., Lomovstev, A., Orientation of collagen at the osteocyte lacunae in human secondary osteons., Journal of Biomechanics, 2008; 41: 3426–3435. doi: 10.1016/j.jbiomech.2008.09.010. van Hove, R.P., Nolte, P.A., Vatsa, A., Semiens, C.M., Salmon, P.L., Smit, T.H., Klein-Nulend, J., Osteocyte morphology in human tibiae of different bone pathologies with different bone mineral density — Is there a role for mechanosensing?, Bone, 2009; 45: 321–329. doi: 10.1016/j.bone.2009.04.238. Britz, H.M., Carter, Y., Jokihaara, J., Leppanen, O.V., Jarvinen, T.L.N., Belev, G., Cooper, D.M.L.,Prolonged unloading in growing rats reduces corticalosteocyte lacunar density and volume in the distal tibia., Bone, 2012; 51: 913–919. doi:10.1016/j.bone.2012.08.112. Mader, K.S., Schneider, Philipp, S., Muller, R., Stampanoni, M.,A quantitative framework for the 3D characterization of the osteocyte lacunar system., Bone, 2013; 57: 142–154. doi: https://doi.org/10.1016/j.bone.2013.06.026. Lanaro, M., Mclaughlin, M.P., Simpson, M.J., Buenzli, P.R., Wong, C.S., Allenby, M.C., Woodruff, M.A., A quantitative analysis of cell bridging kinetics on a scaf- fold using computer vision algorithms., Acta Biomaterialia, 2021; 136: 429–440. doi: https://doi.org/10.1016/j.actbio.2021.09.042. Ruffoni, D., Fratzl, P., Roschger, P., Klaushofer, K., Weinkamer, R.,The bone mineralization density distribution as a fingerprint of the mineralization process., Bone, 2007; 40: 1308– 1319. doi:10.1016/j.bone.2007.01.012. Buenzli, P.R., Osteocytes as a record of bone formation dynamics., Journal of Theoretical Biology, 2015; 364: 418–427. doi:10.1016/j.jtbi.2014.09.028. 20 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Cox, B.N., Purohit, P.K., White, S.N.Strain fields and solitary strain waves as determining factors for the cross-sectional geometry of mouse incisor enamel., Journal of the Mechanics and Physics of Solids, 2015; 193: 105840. doi:10.1016/j.jmps.2024.105840. Guyot, Y., Papantoniou, I., Chai, Y.C., Van Bael, S., Schrooten, J., Geris, L.,A computational model for cell/ECM growth on 3D surfaces using the level set method: a bone tissue engi- neering case study., Biomechanics and Modeling in Mechanobiology, 2014; 13: 1361–1371. doi: 10.1007/s10237-014-0577-5. Callens, S.J.P., Fan, D., van Hengel, I.A.J., Minneboo, M., Díaz-Payno, P.J., Stevens, M.M., Fratila-Apachitei, L.E., Zadpoor, A.A.,Emergent collective organization of bone cells in complex curvature fields., Nature Communications, 2023; 14: 855. doi:10.1038/s41467- 023-36436-w. Ambrosi, D., Guana, F.,Stress-Modulated Growth., Mathematics and Mechanics of Solids, 2007; 12: 319–342. doi:10.1177/1081286505059739. Dunlop,J.W.C.,Fischer,F.D.,Gamsjäger,E.,Fratzl,P., Atheoreticalmodelfortissuegrowthin confined geometries., Journal of the Mechanics and Physics of Solids, 2010; 58: 1073–1087. doi: 10.1016/j.jmps.2010.04.008. Goriely, A.,The Mathematics and Mechanics of Biological Growth., Springer New York, NY, 2017. Joly, P., Duna, G.N., Schone, M., Welzel, P.B., Fruedenberg, U., Werner, C., Petersen, A., Geometry-Driven Cell Organization Determines Tissue Growths in Scaffold Pores: Conse- quences for Fibronectin Organization., PLOS One, 2013; 8: e73545. doi:10.1371/jour- nal.pone.0073545. Gahffarizadeh,A.,Heiland,R.,Friedman,S.H.,Mumenthaler,S.M.,Macklin,P., PhysiCell: An opensourcephysics-basedcellsimulatorfor3-Dmulticellularsystems.,PLOSComputational Biology, 2018; 14: e1005991. doi:10.1371/journal.pcbi.1005991. Alias, M.A., Buenzli, P.R., Modeling the Effect of Curvature on the Collective Behav- ior of Cells Growing New Tissue., Biophysical Journal, 2017; 112: 193–204. doi: 10.1016/j.bpj.2016.11.3203. Hegarty-Cremer,S.G.D.,Simpson,M.J.,Andersen,T.L.,Buenzli,P.R., Modellingcellguidance and curvature control in evolving biological tissues., Journal of Theoretical Biology, 2021; 520: 110658. doi:10.1016/j.jtbi.2021.110658. Kuba, S., Simpson, M.J., Buenzli, P.R., A mathematical model of curvature con- trolled tissue growth incorporating mechanical cell interactions, 2026; BioRxiv. doi: 10.64898/2026.03.10.710423. 21 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Baker, R.E., Parker, A., Simpson, M.J.,A free boundary model of epithelial dynamics., Journal of Theoretical Biology, 2019; 481: 61–74. doi:10.1016/j.jtbi.2018.12.025. Murphy, R.J., Buenzli, P.R., Baker, R.E., Simpson, M.J., A one-dimensional individual- based mechanical model of cell movement in heterogeneous tissues and its coarse- grained approximation., Proceedings of the Royal Society A, 2019; 475: 20180838. doi: 10.1098/rspa.2018.0838. Buenzli, P.R., Kuba, S., Murphy, R.J., Simpson, M.J.,Mechanical cell interactions on curved interfaces., Bulletin of Mathematical Biology, 2025; 87: 29. doi: 10.1007/s11538-024- 01406-w. Brown,P.J.,Green,E.F.,Binder,B.J.,Osborne,J.M., Competingmechanismsforthebucklingof anepithelialmonolayeridentifiedusingmulticellularsimulation.,MathematicalBiosciences, 2025; 380: 109367. doi:https://doi.org/10.1016/j.mbs.2024.109367. Murray, P.J., Edwards, C.M., Tindall, M.J., Maini, P.K.,From a discrete to a continuum model of cell dynamics in one dimension., Physical Review E, 2009; 80: 031912. doi: 10.1103/PhysRevE.80.031912. Murray, P.J., Edwards, C.M., Tindall, M.J., Maini, P.K.,Classifying general nonlinear force laws in cell-based models via the continuum limit., Physical Review E, 2012; 85: 021921. doi: 10.1103/PhysRevE.85.021921. Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., Rubin, D.B.,Bayesian Data Analysis, Third Edition, Chapman and Hall/CRC, 2013. Cleveland, W.S., Grosse, E.,Computational methods for local regression, Statistics and Com- puting, 1991; 1(1), 47-62. doi:10.1007/BF01890836. Purcell, E.J.,Life at low Reynolds number, American Journal of Physics 1977; 45: 3–11. doi: 10.1142/9789814434973_0004. Murphy, R.J., Buenzli, P.R., Baker, R.E., Simpson, S.J., Mechanical Cell Competition in Heterogeneous Epithelial Tissues., Bulletin of Mathematical Biology, 2020; 82: 130. doi: 10.1007/s11538-020-00807-x. Buenzli, P.R., Lanaro, M., Wong, C.S., McLaughlin, M.P., Allenby, M.C., Woodruff, M.A., Simpson, M.J., Cell proliferation and migration explain pore bridging dynamics in 3D printed scaffolds of different pore size., Acta Biomaterialia, 2020; 114: 285–295. doi: 10.1016/j.actbio.2020.07.010. Browning, A.P., Maclaren, O.J., Buenzli, P.R., Lanaro, M., Allenby, M.C., Woodruff, M.A., Simpson, M.J., Model-based data analysis of tissue growth in thin 3D 22 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint printed scaffolds., Journal of Theoretical Biology, 2021; 528: 110852. doi: https://doi.org/10.1016/j.jtbi.2021.110852. Parfitt, A.M.,The physiological and clinical significance of bone histomorphometric data., In Recker R.R. (Ed.), Bone histomorphometry: Techniques and interpretation., 1983. Alias, M.A., Buenzli, P.R.,Osteoblasts infill irregular pores under curvature and porosity controls., Biomechanics and Modeling in Mechanobiology, 2018; 17: 1357–1371. doi: 10.1007/s10237-018-1031-x. Metz, L.N., Martin, R.B., Turner, A.S.,Histomorphometric analysis of the effects of osteo- cyte density on osteonal morphology and remodeling., Bone, 2003; 33: 753–759. doi: 10.1016/S8756-3282(03)00245-X. Hannah, K.M., Thomas, C.D.L., Clement, J.G., Peele, A.G.,Bimodal distribution of osteocyte lacunarsizeinthehumanfemoralcortexasrevealedbymicro-CT.,Bone,2010;47: 866–871. doi: https://doi.org/10.1016/j.bone.2010.07.025. Power, J., Doube, M., van Bezooijen, R.L., Loveridge, N., Reeve, J.,Osteocyte recruitment declines as the osteon fills in: Interacting effects of osteocytic sclerostin and previous hip fracture on the size of cortical canals in the femoral neck., Bone, 2012; 50: 1107–1114. doi: https://doi.org/10.1016/j.bone.2012.01.016. Jones S.P.,Secretory territories and rate of matrix production of osteoblasts., Calcified Tissue Research, 1974; 14: 309–315. doi:10.1007/BF02060305. Martin, R.B., Burr, D.B., Sharkey, N.A.,Skeletal Tissue Mechanics., Springer New York, NY, 1998. Buenzli, P.R., Pivonka, P., Smith, D.W.,Bone refilling in cortical basic multicellular units: insights into tetracycline double labelling from a computational model., Biomechanics and Modeling in Mechanobiology, 2014; 13: 185–204. doi:10.1007/s10237-013-0495-y. Bidan, C.M., Kollmannsberger, P., Gering, V., Ehrig, S., Joly, P., Petersen, A., Vogel, V., Fratzl, P., Dunlop, J.W.C.,Gradual conversion of cellular stress patterns into pre-stressed matrix architecture during in vitro tissue growth., Journal of The Royal Society Interface, 2016; 13: 20160136. doi:10.1098/rsif.2016.0136. VandenHeuvel,D.J.,Devlin,B.L.,Buenzli,P.R.,Woodruff,M.A.,Simpson,M.J., Newcomputa- tionaltoolsandexperimentsrevealhowgeometryaffectstissuegrowthin3Dprintedscaffolds., Chemical Engineering Journal, 2023; 475: 145776. doi:10.1016/j.cej.2023.145776. Lassen, N.E., Andersen, T.L., Pløen G.G., Søe, K., Hauge, E.M., Harving, S., Eschen, G.E.T., Delaisse, J-M.,Coupling of Bone Resorption and Formation in Real Time: New Knowledge 23 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint Gained From Human Haversian BMUs., Journal of Bone and Mineral Research, 2017; 32: 1395–1405. doi: 10.1002/jbmr.3091. Robling, A., Stout, S.,Morphology of the drifting osteon., Cells Tissues Organs, 1998; 164: 192–204. doi:10.1159/000016659. Schnitzler, C.M., Mesquita, J.M.,Cortical porosity in children is determined by age-dependent osteonal morphology., Bone, 2013; 55: 476–486. doi:10.1016/j.bone.2013.03.021. Cooke, K.M., Mahoney, P., Miszkiewicz, J.J.,Secondary osteon variants and remodeling in human bone., The Anatomical Record, 2022; 305: 1299–1315. doi:10.1002/ar.24646. Hegarty-cremer, S.G.D., Borggaard, X.G., Andreasen, C.M., van der Eerden, B.C.J., Simp- son, M.J., Andersen, T.L., Buenzli, P.R.,How osteons form: A quantitative hypothesis- testing analysis of cortical pore filling and wall asymmetry., Bone, 2024; 180: 116998. doi: 10.1016/j.bone.2023.116998. Tee, S.Y., Bausch, A.R., Janmey, P.A.,The mechanical cell., Current Biology, 2009; 19(17): R745 - R748. doi:10.1016/j.cub.2009.06.034. Marotti, G., Zallone, A.Z., Ledda, M., Number, Size and Arrangement of Osteoblasts in Osteons at Different Stages of Formation., Calcified Tissues, 1976; 21: 96–101. doi: 10.1007/BF02546434. Werner, M., Blanquer, S.B.G., Haimi, S.P., Korus, G., Dunlop, J.W.C., Duda, G.N., Grijpma, D.W., Petersen, A., Surface Curvature Differentially Regulates Stem Cell Migration and Differentiation via Altered Attachment Morphology and Nuclear Deformation., Advanced science, 2017; 1600347. doi:10.1002/advs.201600347. Gudipaty, S.A., Lindbolm, J., loftus, P.D., Redd, M.J., Edes, K., Davey, C.F., Krishnegowda, V., Rosenblatt, J.,Mechanical stretch triggers rapid epithelial cell division through piezol., Nature, 2017; 118-121. doi:10.1038/nature21407. Ravasio, A., Cheddadi, I., Chen, T., Pereira, T., Ong, H.T., Bertocchi, C., Brugues, A., Jacinto, A., Kabla, A.J., Toyama, Y., Trepat, X., Gov, N., Neves de Almeida, L., Ladoux, B.,Gap geometry dictates epithelial closure efficiency., Nature communications, 2015; 6, 7683. doi: 10.1038/ncomms8683. Kuba, S., Simpson, M. J., and Buenzli, P. R. (2026). Stochastic model of curvature- controlled tissue growth incorporating mechanical cell interactions. GitHub repository. https://github.com/Shahak-Kuba/2026_Stochastic_Tissue_Growth 24 .CC-BY 4.0 International licenseperpetuity. It is made available under a preprint (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in The copyright holder for thisthis version posted March 16, 2026. ; https://doi.org/10.64898/2026.03.12.711199doi: bioRxiv preprint