Understanding how tissue growth in porous sca ffolds is influenced by geometry is a fundamental challenge in the field of tissue engineering. We investigate the influence of pore geometry on tissue growth using osteoblastic cells in 3D printed melt electrowritten sca ffolds with square-shaped pores and non-square pores with wave-shaped boundaries. Using a reaction-di ffusion model together with a likelihood-based uncertainty quantification framework, we quantify how the cellular mechanisms of cell migration and cell proliferation drive tissue growth for each pore geometry. Our results show that the rates of cell migration and cell proliferation appear to be largely independent of the pore geometries considered, suggesting that observed curvature e ffects on local rates of tissue growth are due to space availability rather than directly a ffecting cell behaviour. This result allows for simple squared-shaped pores to be used for estimating parameters and making predictions about tissue growth in more realistic pores with more realistic, complicated shapes. Our findings have important implications for the development of predictive tools for tissue engineering and experimental design, highlighting new avenues for future research. Keywords: Tissue engineering; Mathematical modelling; Curvature; Tissue growth; Pore geometry 1. Introduction Tissue engineering aims to regenerate damaged or diseased tissues [1]. A key challenge in tissue engineering is to understand how various clinically-motivated experimental conditions influence tissue growth [2]. Recent ad- vancements in three-dimensional (3D) printing allow us to investigate tissue growth in 3D-printed sca ffolds of various shapes and sizes, enabling realistic migration and proliferation behaviours to be studied in well-controlled experi- mental conditions [3, 4]. While the e ffect of pore geometry and tissue curvature on tissue growth is well-known [5–7], understanding how these e ffects relate to cellular-level mechanisms remains poorly understood. Understanding these cellular mechanisms would enable the prediction and analysis of tissue growth in complex geometries from the calibration of mathematical models in simpler geometries, providing a valuable computational tool for screening experimental designs [8] of sca ffold geometries and providing plausible results on these new geometries. This kind ∗Corresponding author. Email: matthew.simpson@qut.edu.au Preprint submitted to Elsevier August 21, 2023 of computational tool could enable more personalised approaches to tissue engineering, in which a specific sca ffold size and shape could be tailored to an individual patient, making use of these predictions to screen for possible patient outcomes. 3D printing technology for biofabrication has evolved rapidly [9, 10]. Melt electrowriting, a technique for high quality 3D printing, allows for micro- and nano-scale fibres to created, enabling great control over the fibres and the pore geometry, making it possible to produce realistic sca ffold geometries with a regular array of pores for growing tissue [11–14]. These sca ffolds are designed so that cells and tissues experience a similar mechanical support as they would, for example, in skin and bone tissues, enabling realistic cell migration and cell proliferation behaviours to be observed and measured [3]. Previous work has focused primarily on squared-shaped pores [15, 16], although more complicated sca ffolds can also be produced [10, 17]. One important factor in understanding tissue growth is curvature [18, 19]. In the context of bone tissue, Bidan et al. [20] suggest that cell tension can influence tissue curvature that, in turn, can stimulate tissue growth. Callens et al. [21] discuss how cells respond to their surrounding geometry, even across large spatial scales, and how this a ffects bone tissue growth. Mathematical modelling studies performed by Alias and Buenzli [22–24] and Hegarty-Cremer et al. [25] also investigate the role of geometry and curvature on tissue growth and cell crowding in bone tissue growth experiments. Figure 1: Sca ffold pore boundary (red boundary), void (black interior), void boundary (magenta boundary), tissue (blue /green region), and fibres (exterior black boundary). The white outlines show the sca ffold boundaries. The blue channel in the microscope images shows the cell nuclei (DAPI), and the green channel shows the tissue and cytoskeleton (phalloidin). The DAPI on the right image is shown in grey. Tissue growth experiments in porous sca ffolds are of great importance in tissue engineering [16, 26]. In these experiments, cells are seeded onto the perimeter of a sca ffold, leading to cell migration and cell proliferation that produces an inward-growing tissue. The shape of the region that is devoid of cell and tissue material, referred to as thevoid in Figure 1, matches the shape of the sca ffold boundary for early times, rounding o ffover time, eventually 2 forming a circular front until the void closes, which we call pore bridging . The feature of interest in these experi- ments is the time that the tissue bridges, called the bridging time [16, 26]. A circular front arises with time in many pore sca ffold geometries such as square, triangular, and hexagonally shaped pores [15, 16, 27, 28], though the precise mechanisms driving the cells into these circular fronts in general geometries remains unclear. The most common way to report a set of pore bridging experiments is to record snapshot images of the tissue growth process [15, 16, 26]. This approach allows us to estimate the bridging time within an interval instead of identifying the precise time of bridg- ing [16, 26], thereby introducing some uncertainty into the experimental estimate of the bridging time. To interpret such measurements meaningfully, we are interested in developing mathematical modelling tools that mechanistically capture cell migration and cell proliferation within a framework that explicitly incorporates uncertainty in the ex- perimental measurements, uncertainty in the parameter estimates in the mathematical model, and that is capable of making predictions of new experiments that incorporate these uncertainties. This will allow us, for example, to study how variability in bridging times can be explicitly integrated into the mathematical model, as well as interpreting the predictions of the mathematical model. In this paper, we study mechanistic cell behaviour in pore bridging experiments performed within pores on 3D- printed sca ffolds made from polycaprolactone [4, 10, 26]. We consider tissue growth in two di fferent shaped pores, a square-shaped pore (Figure 1(a)) and a wave-shaped pore (Figure 1(b)), with the aim to understand whether the details of cell migration and cell proliferation are a ffected by di fferences in the pore geometry. In particular, we are interested in determining whether the cellular mechanisms driving tissue growth in the more realistic wave-shaped pore are indistinguishable from the cellular mechanisms driving tissue growth in the simpler square geometry. Moreover, we are interested in obtaining estimates of bridging time with uncertainty, for both pore geometries, through probability distributions that allow a user to predict probabilities of bridging time occurring within a specified time window for any pore geometry. All experiments reported in this paper involve tissue growth using murine calvarial osteoblastic cells (MC3T3-E1) [29]. Our group has previously used these experiments on the square pores with di fferent sizes [15, 16, 26], investigating the relationship between the cell migration rate measured in terms of the cell di ffusivity D, cell proliferation rate λ, and sca ffold size. This previous work showed that the product Dλ, which controls the long-time rate of tissue production, appeared to be una ffected by the pore size, but did not consider the role of pore shape. These previous experiments displayed variability in the time to bridge, even in well controlled experiments with the same pore size. This variability motivates the need for developing mathematical modelling tools that incorporate variability and uncertainty quantification into predictions. Hence, our analysis uses a combination of numerical simulations from a mathematical model and with statistical analysis that takes numerical simulations and quantifies results together with uncertainty. These considerations will be used to answer the following broad questions: 1. Are the cellular mechanisms driving tissue growth independent of pore shape? 2. Can we use results on one pore shape to make predictions, with uncertainty, on another geometry? 3. What data, and how much data, is su fficient for accurately comparing results for di fferent shaped pores? 3 We address these questions using a model-based approach, using the Porous-Fisher partial di fferential equation (PDE) as a model for tissue growth driven by combined cell migration and cell proliferation [30]. Solving this mathematical model requires parameter estimates for the cell di ffusivity, D, the cell proliferation rate, λ, and the initial cell density on the sca ffold boundary, u0. Estimates for Dandλcannot be obtained directly from experimental images, hence methods that use information about the images are needed. For each pore shape, we calibrate this mathematical model with an experimental dataset containing information about the position of the tissue front over time. We apply a likelihood-based analysis to this dataset [31] with the aim to estimate the combined e ffect of proliferation and migration rates, by estimating the product Dλ. Profile likelihoods are used to quantify the uncertainty in Dλ, providing confidence intervals for Dλ[32], and allowing us to determine what parameters or parameter combinations can be estimated [33, 34], providing insights into the second and third research questions listed above. The confidence intervals obtained on each geometry can be used to compare the tissue growth mechanisms for each pore geometry to answer the first research question. This analysis also enables us to make predictions with uncertainty about the pore bridging time. By using this likelihood-based approach to make predictions, we can take results from the square- shaped pores and estimate, with uncertainty, the pore bridging times on the wave-shaped pores. Similarly, we can use the wave-shaped pores to make predictions on the square-shaped pores. Comparing both situations allows us to answer the second research question. The simplest interpretation of our results is that the cell migration and cell proliferation rates are independent of the pore sca ffold geometry. 2. Materials and Methods In this section, we describe the methods used for the experiments and the data that we collect from these ex- periments. Following this description, we introduce the mathematical models we use and how we apply likelihood analysis for performing statistical inference from these experiments. 2.1. Tissue growth experiments All reagents are sourced from Thermo Fisher unless otherwise stated. Using a melt electrowriting printer described previously [26], polycaprolactone (45 kDa, Sigma Aldrich) fibres of diameter 50 µm are fabricated into a three-layer scaffold which was then biopsy punched to form a 6 mm disc. The resultant sca ffold has a thickness of approximately 150µm. The code used to produce the outline of each pore is adjusted to produce square-shaped and wave-shaped pores of comparable size, both being derived from a unit of cell of 500 µm. Prior to cell seeding, sca ffolds are sterilised under UV light overnight. The cells used are murine calvarial osteoblastic cells (MC3T3-E1) [29] that are cultured in α-MEM, 10 % fetal bovine serum, and 1 % penicillin-streptomycin, and are approximately 20–30 µm in diameter [26]. Cells were expanded in a T75 culture flask and at 80 % confluency were detached with TrypLE. Cells were seeded at a density of 10,000 cells per sca ffold in 48-well plates. After allowing 4 h for the cells to attach to each sca ffold after seeding, an additional 500 µL was added. Cell-seeded sca ffolds were cultured in a humidified environment at 37◦C in 5 % CO2 for 28 days. The media is changed every 2–3 days from day 5 to day 14, every 1–2 4 days from day 15 to 28. The viability of the cells are assessed at day 10, 14, and 28 using calcein-AM and ethidium homodimer-1 to stain live and dead cells, respectively. At specific timepoints, cell-seeded sca ffolds were fixed with 4 % paraformaldehyde and stained with 4’,6-diamidino-2-phenylindole (DAPI) and Alexa Fluor ™488 Phalloidin, which stain cell nuclei and actin filaments, respectively. High resolution images of the centre of each sca ffold are obtained using fluorescent microscopy (Zeiss, AxioObserver 7). For each pore shape and timepoint, fixation, staining, and microscopy are repeated across two or three identically prepared replicates. Each experimental replicate provides information of several pores, giving information about tissue growth data from day 5 to day 28. This procedure gives 41 data points for the square pore, and 3 data points for the wave pore; while this number of data points is relatively small, it is su fficient for our analysis as we use a likelihood-based uncertainty quantification framework that naturally incorporates variability associated with finite sample sizes, as introduced in Section 2.4. More information on the procedure for capturing these images is provided in [26]. These experiments lead to a tissue composed of a monolayer of cells within the pore. The growth of monolayers is achieved through the design of the pore geometry as the width of the pores is large compared to the pore depth. The vertical pore depth is approximately equal to the average cell diameter, 20–30 µm [26], which means that we do not observe cells growing on top of each other in the vertical direction. 2.2. Data and image processing The experiments provide us with several images at four time points (days 7, 14, 25, and 28), and each image contains information about several pores. To summarise the tissue growth processes, for each pore we calculate two quantities: yi,j c=area of void area of pore,yi,j p=perimeter of void perimeter of pore, (1) where yi,j candyi,j pdenote the void coverage andnormalised void perimeter , respectively, for the jth pore at the ith time, ti. We do not collect data from any pores that have bridged, or pores that have bridged in a way so that the void splits into multiple disconnected regions. All images at t=25 day and t=28 day show that all pores have bridged before 25 days, so we focus on measuring (1) at t1=7 day and t2=14 day. A precise description of how we compute the quantities in (1) from the images is given in Section S1 of the Supplementary Material. To complete the processing, we select the computational representations for the square and wave geometries for use in the mathematical model described Section 2.3 (Figure 2) for comparison with each data point. In the square case, we construct this independent of the images, and simply define a boundary for a square with its lower-left corner at the origin and side length L=475µm. For the wave geometry, we take a single image from the experiments and choose its boundary as a representative boundary for each experiment, which is reasonable as all sca ffolds are uniformly printed. 5 2.3. Mathematical model We use the Porous-Fisher PDE to model tissue growth, as this model explicitly describes how cell migration and cell proliferation leads to tissue growth with sharp fronts that we observe in the experiments (Figures 3 and 4) [15, 30]. Since the pore bridging process starts after 5 days, the PDE is solved for t>5 day. Letting ˜ u(x,y,t) [cells /µm2] denote the density of cells at a point ( x,y) and time t[day], and ˜Kdenoting the maximum carrying capacity density [cells /µm2], we define a normalised density u(x,y,t)∈[0,1] by u(x,y,t)=˜u(x,y,t)/˜K. Thus, the model for u(x,y,t) is given by ∂u(x,y,t) ∂t=D(∂ ∂x" u(x,y,t)∂u(x,y,t) ∂x# +∂ ∂y" u(x,y,t)∂u(x,y,t) ∂y#) | {z } contact stimulated cell migration+λu(x,y,t)1−u(x,y,t) | {z } contact inhibited cell proliferation,(x,y)∈Ω,(2) du(x,y,t) dt=λu(x,y,t)1−u(x,y,t) | {z } contact inhibited cell proliferation,(x,y)∈∂Ω, (3) u(x,y,5)=u0(x,y)∈∂Ω, 0 ( x,y)∈Ω.(4) where (2) is applied on the interior sca ffold pore space Ω, the space inside the red curves of Figure 2; (3) is applied on the boundary ∂Ω, the red curve in Figure 2. We note that while our sca ffolds are three-dimensional, their thickness is small compared to their width, and so this two-dimensional depth-averaged model (2)–(4) is reasonable [35]. This model is characterised by three parameters ( D,λ,u0), where D[µm2/day] is the cell di ffusivity that controls the rate of cell migration, λ[day−1] is the cell proliferation rate, and u0is the normalised density of cells on the sca ffold boundary ∂Ωatt=5 day. We solve Equations (2)–(4) numerically using the finite volume method with an unstructured triangular mesh (Figure 2), as described in Section S2 of the Supplementary Material. Figure 2: Schematics for triangular meshes used for numerical solutions of Equations (2)–(4) on the (a) square-shaped pore, (b) the wave-like pore, and (c) the cross-shaped pore considered in Section 3.3. Denser meshes are used for the actual solutions. The triangles represent the mesh we use for studying these geometries computationally, as discussed in Section 2.3. The red curves represent the boundary ∂Ω, and the region bounded by these curves is the interior domain Ω. 6 Solutions of Equations (2)–(4) are used to compute model predictions corresponding to the data yi,j candyi,j pfrom (1). Using the numerical solution for u(x,y,t), we identify the contour u(x,y,t)=1/2 to indicate the location of the tissue front [15, 16]. We define predictions of the yi,j candyi,j pby µc(t)=1 A(Ω)nX k=1Ak(t), µ p(t)=1 ℓ(∂Ω)nX k=1ℓk(t), (5) respectively, where nis the number of triangular elements in the mesh, Ak(t) is the area of the portion of the kth element at the time tthat is inside the contour u(x,y,t)=1/2,ℓk(t) is the length of the line through the kth element at the time tthat is on the contour u(x,y,t)=1/2 or zero if the contour does not go through the element, and A(Ω) and ℓ(∂Ω) are the area and perimeter of the domain Ω, respectively. More details on how we compute the coverage µcand normalised perimeter µpare given in Section S3 of the Supplementary Material. We remark that our definition of the PDE (2)–(4) involves working with a nondimensional dependent variable u and nondimensional parameter u0, while retaining the dimensional parameters Dandλand dimensional variables x,y, andt, so that spatial and temporal features can be compared with experimental images, similar to [36]. Our interpreta- tion of the dependent variable, u(x,y,t), is di fferent, though. While it is possible to work with the dimensional density ˜u(x,y,t) [15, 37], this would require manually counting cells to estimate cell densities in space and time [36, 38, 39]. Thus, to be consistent with the fact that we only treat the leading edge in the experimental images for computing yi,j c andyi,j p, it makes sense to consider the dimensionless ratio u(x,y,t)=˜u(x,y,t)/˜Krather than ˜ u(x,y,t) itself. 2.4. Parameter estimation We use a likelihood-based approach to estimate model parameters [33, 40]. This approach takes predictions µc(t) andµp(t) from solutions of Equations (2)–(4) and compares them with the noisy experimental observations, yi,j cand yi,j p. We assume that these noisy experimental observations are all independent realisations of random variables Yi cand Yi p, respectively, that are defined by [15] Yi c∼N µc(ti;θ),σ2 c ,and Yi p∼N µp(ti;θ),σ2 p , (6) whereN(µ,σ2) denotes the normal distribution with mean µand variance σ2. These random variables in (6) each have a mean that depends on θ=(D,λ,u0), and the variances σ2 candσ2 pneed to be estimated. The values of σc andσpare treated as constants that are pre-estimated using the sample standard deviation of the experimental data aggregated for each tiand each j, as described in Section S3 of the Supplementary Material. Given su fficient experimental data we could, in theory, estimate all model parameters θandσ2 candσ2 pdirectly, however, as we will show, our data is insu fficient for this purpose. We find that our numerical simulations of Equations (2)–(4) on the time scale of our experiments are relatively independent of u0, and so we show results for a range of pre-specified values of u0rather than focusing on any single value, demonstrating this independence. While it would be ideal to estimate both Dandλ, we find that it is di fficult to treat them separately, as we show in Section S4 of the Supplementary Material [16, 37]. For our purposes, though, we are mainly interested in the combined e ffectDλ, 7 as this is the variable that a ffects the velocity of the tissue and thus the bridging time [16]. Thus, rather than using θ=(D,λ,u0), we instead re-parametrise the vector of model parameters as θ=(Dλ,λ), omitting u0and following [41]. With this definition, µc(ti;θ) andµp(ti;θ) still refer to predictions from the solution of the PDE with parameters Dandλ, and the choice of the fixed value of u0is left implicit. 2.4.1. Log-likelihood function The log-likelihood function ℓ(θ|y) is a function of the model parameters that describes the likelihood that the model has parameter values θgiven that the data observed is y[31, 42]. In this work we have ℓ(θ|y)=2X i=1J(i)X j=1h logϕ yi,j c;µc(ti;θ),σ2 c +logϕ yi,j p;µp(ti;θ),σ2 pi , (7) where yis the vector of observations, ϕ(x;µ,σ2)=(2πσ2)−1/2exph −(x−µ)2/(2σ2)i is the normal probability density function, and J(i) is the number of pores included at the time ti. For the wave pore, the log-likelihood (7) only includes the sum at i=2 as there is only data at t2=14 day in this case. For notational convenience we write ℓ(θ|y) as ℓ(θ). The log-likelihood depends on the parameters θ=(Dλ,λ) and the fixed values for u0,σc, andσp, with{D,λ,u0} governing the solution to the PDE and σpandσcgoverning the measurement model. 2.4.2. Maximum likelihood estimation We obtain a best-fit estimate for the parameters θby maximising ℓ(θ). This procedure is called maximum likeli- hood estimation [31], and it results in a maximum likelihood estimate (MLE) for θ, denoted ˆθ. For this maximisation, we constrain the values of Dandλso that 0 µm2/day2<Dλ⩽10 000 µm2/day2. Forλwe use 0 day−1<λ⩽5 day−1 in the square and 0 day−1< λ⩽10 day−1in the wave. These bounds do not a ffect the results significantly. More detail on how we perform this maximum likelihood estimation is given in Section S5 of the Supplementary Material. 2.4.3. Uncertainty quantification One limitation of maximum likelihood estimation is that we obtain a single point estimate for the MLE ˆθ, and the asymptotic uncertainty in this point estimate depends upon the curvature of the log-likelihood function [43]. To quantify the uncertainty in this estimate, we combine two approaches. For the first approach, we evaluate the log-likelihood function ℓ(θ) over a large grid of Dλandλvalues. We then use this grid to find all points where ℓ(θ)−ℓ∗⩾−χ2 2,1−α/2, whereχ2 d,qis the qth quantile of the χ2distribution with ddegrees of freedom and ℓ∗=ℓ(ˆθ) is the maximum likelihood, as the resulting set of values defines a 100(1 −α)% credible region (CR) for θ[31, 44]. We useα=0.05 in this work, giving −χ2 2,0.95/2≈−3. The approach above gives us a two-dimensional region representing the uncertainty in θ. It will also be useful to reduce these regions to confidence intervals for each parameter, and most importantly for the parameter combination Dλ, which allows us to make predictions about the variability in tissue growth later in Section 2.4.5. We take a profile likelihood approach to consider each parameter individually, specifying a range of values for an interest parameter 8 and using numerical optimisation, reducing the log-likelihood function to a series of interpretable univariate func- tions. These univariate results then provide insight into the curvature of the log-likelihood function, and hence the uncertainty in the MLE point estimates. The resulting univariate function results in what is known as the profile log- likelihood function [31]. Following [15, 33, 40], we define the profile log-likelihood in terms of the interest parameter Dλ. For a given value of Dλ, we define the normalised profile log-likelihood: ℓp(Dλ)=max λ∈Λ[ℓ(Dλ,λ)]−ℓ∗, (8) whereℓ(Dλ,λ)=ℓ(θ),ℓ∗=ℓ(ˆθ), and Λ ={λ: 0 day−1< λ⩽5 day−1}. This definition gives a simple univariate function of Dλthat reaches a maximum of zero at the MLE, and the curvature of this function is related to inferential precision — a profile log-likelihood function with a well-defined peak at zero indicates a parameter that has been well estimated and identified, while a flat profile means that the data was insu fficient for estimating or obtaining any inference for that parameter [40, 45]. A useful feature of (8) is that it can be used for constructing approximate confidence intervals for Dλ, with an approximate 100(1 −α)% confidence interval given by the set of all Dλsuch that ℓp(Dλ)⩾−χ2 1,1−α/2 [31]. In this work, we use α=0.05 so that we are constructing 95% confidence intervals, giving c∗=−χ2 1,0.95/2≈−1.92. The procedure we use for computing profile likelihoods is implemented in the J ulia package ProfileLikelihood.jl [46], and a summary of the procedure is outlined in Section S5 of the Supplementary Material. 2.4.4. Parameter-wise prediction intervals The two-dimensional likelihood function allows us to propagate the uncertainty in Dλthrough to give us a pre- diction interval in terms of the outcome of the mathematical model for a variable of interest, such as the bridging time or cell density. Using the approach developed by [33, 40] which builds on basic properties of likelihood function [31], we are able to quantify the uncertainty in the cell densities u(x,y,t) and the bridging time tbdirectly from our likelihood function. In particular, by taking pairs of parameter values inside of the 95% credible region from the log-likelihood function and computing the variable of interest at each pair, we obtain a sample of values that gives the uncertainty in our variable of interest, as described in Section S5 of the Supplementary Material. We use this method to obtain prediction intervals for yi,j candyi,j pover time. Moreover, we can obtain prediction intervals for the bridging time, tb, at whichµc(t) first becomes zero. This procedure returns, in addition to the prediction intervals, a sample of bridging times, which we use to obtain probability distributions for the bridging time via KernelDensity.jl [47]. We represent this probability distribution using a probability density function (PDF) p(tb) for the bridging time which can be understood as [42] p(tb)∆t≈P(t<tb<t+ ∆t), (9) where P(t<tb<t+ ∆t) is the probability that the bridging time tbis between tandt+ ∆t, given the uncertainty in the parameters θ, and ∆tis some small su fficiently interval of time. This PDF p(tb) allows us to compute probabilities that the bridging time occurs in anygiven interval. The complete procedure for how we obtain these results is implemented 9 inProfileLikelihood.jl [46], and the method that we implement is outlined in Section S5 of the Supplementary Material. 2.4.5. Predicting variability in tissue growth We further extend our results to provide a more qualitative approach to assess the uncertainty in the tissue growth on these pores. Taking values for the parameters inside their confidence intervals from the profile likelihoods, we can produce time series model predictions of the solution to Equations (2)–(4), indicating the variability that we might expect in the tissue growth. For these predictions, we take three parameter values for θ: (1) ˆθ, the MLE; (2) ˆθL, where we take Dλto be the lower endpoint of its confidence interval from ℓp(Dλ) andλto be the lower endpoint of its confidence interval from ℓp(λ); (3) ˆθU, where we take Dλto be the upper endpoint of its confidence interval from ℓp(Dλ) andλto be the upper endpoint of its confidence interval from ℓp(λ). Solving Equations (2)–(4) with these three combinations of θprovides a simple way of giving a visual interpretation of the uncertainty in cell density as an approximation to the true uncertainty bounds. 3. Results We now give the results from our experiments and from our likelihood analysis. Following these results, we conclude with a description of how we can use these results to predict future pore bridging experiments. 3.1. Experimental images A subset of the results for the pore bridging experiments on the square geometry are shown in Figure 3 for days 7, 14, 25, and 28, where we see most of the pores take longer than 14 days to bridge, although there is some significant variability in this bridging time as we can even see some pores have completely bridged by day 14. In each pore, the growing tissue always forms a circle before bridging. In total, we have n=41 imaged pores included in the dataset for the square, with 26 at day 7 and 15 at day 14. The results we use for the pore bridging experiments on the wave geometry come only from day 14, and they are shown in Figure 4. Just as we saw in Figure 3, there is significant variability in the bridging time – while most pores appear to have bridged by day 14, some are still open, with a few being far from closed. In these wave pores, the void appears to initially close in as an oval before the void boundary eventually forms a circle. Our interest is in comparing the cell migration and cell proliferation rates between the pores of Figure 3 and Figure 4, but from these images it is not immediately clear whether these are similar or not. Section 3.2 shows results making this comparison using a mathematical model. In total, we only have n=3 pore images for the wave geometry, all at day 14, since all other pores are closed and thus no other data is available for yi,j candyi,j p, or parts of the sca ffolds were not imaged as in the leftmost pore in the first image of Figure 4. The pores used for the data are given in the first, fourth, and sixth images in Figure 4. 10 Figure 3: Experimental images for the square geometry. The images are composite fluoresence microscopy images of pore bridging experiments, with the blue channel showing the cell nuclei (stained with DAPI); the green channel showing the tissue and cytoskeleton (stained with phalloidin). Note that each image is from an independent experiment. Figure 4: Experimental images for the wave geometry on day 14 of each experiment. The images are composite fluoresence microscopy images of pore bridging experiments, with the blue channel showing the cell nuclei (stained with DAPI); the green channel showing the tissue and cytoskeleton (stained with phalloidin). Note that each image is from an independent experiment. 11 3.2. Parameter estimation and parameter identifiability We now consider the likelihoods, profile likelihoods, prediction intervals, and tissue growth predictions for the square and wave pores, demonstrating how well we can calibrate our model to the experimental data and make predictions between the two geometries. 3.2.1. Square pore D𝜆 [𝜇m2/day2]0100200300 400𝜆 [1/day] 0.02.55.0(a): Likelihood CRs D𝜆 [𝜇m2/day2]0100200300 400ℓp(D𝜆) −2−10(b): Profiles for each u0t [day]0 10 20 3 0 40𝜇c(t) 0.00.250.50.751.0(c): Square coveraget [day]0 10 20 30 40𝜇p(t) 0.00.250.50.751.0(e): Square perimetertb [day]22 2 6 3034p(tb) 0.00.1(g): Square p (tb) t [day]0 10 20 3 0 40𝜇c(t) 0.00.250.50.751.0(d): Wave coveraget [day]0 10 20 30 40𝜇p(t) 0.00.250.50.751.0(f): Wave perimetertb [day]19 2 3 2 7p(t b) 0.00.10.2(h): Wave p(tb)u0 0.2 0.3 0.4 Figure 5: Likelihood analysis results for the square pore. In (a), the lines give the boundaries of the 95% credible region for θfor each u0, and the vertical dashed lines show the MLE for Dλ(see Table 1). The profile log-likelihoods for Dλfor each u0are shown in (b), with the threshold c∗≈−1.92 shown with a horizontal red line and the vertical dashed lines show the MLEs for Dλ. In (c)–(f), predictions for µc(t) andµp(t) on each pore geometry are shown, with the blue dots showing the experimental data, the surrounding solid lines giving 95% prediction intervals for each u0, and the dashed lines showing the corresponding estimates at the MLE ˆθ. The estimates for the PDF p(tb) of the bridging time on each pore geometry are shown in (g)–(h). The results in (d), (f), and (h) are predictions on the wave geometry using parameters inferred from the square pore data. Figure 5(a) shows that the credible regions for ( Dλ,λ) for each u0have a similar shape, and the MLEs for Dλare all around the same value. The boundary of each credible region is well-defined in Dλbut not forλ, indicating that we are only able to obtain reliable estimates for Dλbut not for λas we might anticipate [16, 37]. The confidence intervals we obtain from the profile log-likelihoods shown in Figure 5(b) for Dλare approximately the same for each u0, given approximately by 90 µm2/day2<Dλ < 300µm2/day2for each u0. The width of this interval is relatively small, noting that previously reported estimates of Din the literature vary across several orders of magnitude [15, 48]. The predictions for µc(t) on each geometry are shown in Figure 5(c)–(d). We see in (c) that we can recover the data on the square, with the prediction intervals capturing the average experimental data points, and the prediction intervals are indistinguishable for each u0. The dashed lines show the predictions from the MLE ˆθ, indicating the most likely outcome of the experiments, and these curves too pass through the average of the experimental data points and are indistinguishable for each u0. In contrast, we see for µc(t) on the wave geometry that we do not capture the precise values for yi,j c, although if we had more data points then we would likely capture more values due to their variability. These results are also independent for u0. The corresponding figures for µp(t) are shown in Figure 5(e)–(f), where we again capture the data on the square pore but not the data from the wave pore, and again the curves are all independent 12 ofu0. Lastly, we show the probability distributions for the bridging time (9) on each pore geometry in Figure 5(g)– (h). These distributions have a similar shape for each u0. The mode for tbon the square pore appears to be around tb≈24 days, and the distribution shows that we expect more pores to bridge between 23 to 30 days, consistent with our experiments in Figure 3. Similarly, the mode for tbis around 20 days on the wave geometry, with most pores bridging within 19 to 24 days. The observation that the computed results in Figure 5, and in Figures 6–10, are all indistinguishable for each u0is an important result – it suggests that both the mean and variability in our estimates are relatively insensitive to u0, indicating that precise measurements of u0are not critical. y [𝜇m] 0200400Day 7Day 1 4Day 2 5Day 28y [𝜇m] 0200400 x [𝜇m]0 200 4 00y [𝜇m] 0200400 x [𝜇m]0 200 40 0x [𝜇m]0 2 00 400x [𝜇m]0 200 4 00𝜃̂L 𝜃̂ 𝜃̂U Figure 6: Model predictions for the variability in the tissue growth behaviour for the square geometry for u0=0.2 at days 7, 14, 25, and 28. The parameters used are ˆθL=(ˆDL,ˆλL)=(73µm2/day,1.2 day−1),ˆθ=(ˆD,ˆλ)=(32µm2/day,4.9 day−1), and ˆθL=(ˆDU,ˆλU)=(60µm2/day,5 day−1). Figure 6 shows a summary of model predictions where we explore the variability in the behaviour of the experi- ments over time. In particular, we show numerical solutions of Equations (2)–(4) at ˆθ,ˆθL, and ˆθU, each for u0=0.2. In this figure, the columns show predictions for a given time, with each row corresponding to a di fferent value for θ. The middle row corresponds to ˆθ, meaning the prediction that we expect to be most likely. The first two columns show model predictions for the two days that we use for calibrating the model, while the last two columns are genuine 13 predictions since our mathematical model is not calibrated to data from these predictions. We see that there is a lot of variability in the position of the void boundary, especially at day 14, depending on the choice of ˆθL,ˆθ, orˆθU, which approximately matches the observed variability in the experimental images in Figure 3. The numerical results on day 14 show that, at θ=ˆθL, the void boundary is still close to the pore boundary, but the bottom row shows that, at θ=ˆθU, the pore is half-way to being filled. The tissue boundaries do not round o ffas clearly as in Figure 3, although there is some rounding in the corners of these boundaries; this is a limitation of how we calibrate our model, where we are instead seeing results averaged over more shapes. For some simulations in between the limits of the confidence intervals, the tissue boundary rounds o ffmore strongly (not shown). y [𝜇m] 0200400600Day 7Day 1 4Day 2 5Day 28y [𝜇m] 0200400600 x [𝜇m]0200 400600y [𝜇m] 0200400600x [𝜇m]0200400600x [𝜇m]0200 400 600x [𝜇m]0 2 00400600𝜃̂L 𝜃̂ 𝜃̂U Figure 7: Model predictions for the variability in the tissue growth behaviour for the wave geometry using results from the square geometry for u0=0.2 at days 7, 14, 25, and 28. The parameters used are ˆθL=(ˆDL,ˆλL)=(73µm2/day,1.2 day−1),ˆθ=(ˆD,ˆλ)=(32µm2/day,4.9 day−1), and ˆθL=(ˆDU,ˆλU)=(60µm2/day,5 day−1). The model predictions of the tissue growth in Figure 7 show predictions of how the wave pores will evolve over time for u0=0.2, using parameter estimates obtained by calibrating Equations (2)–(4) to data from the square pores. Similar to what we noted in Figure 6, we do not see the same circular voids in Figure 7 as we do in the experimental images in Figure 4, though this is expected as the predictions average over many curves. Similarly, we see high 14 variability in the results, with pores at day 14 ranging from being slightly closed to more than half-way closed. This variability is a positive result, matching the variability in the experimental images (Figure 4). Together, these model predictions indicate that our model has been well-calibrated to the experimental data on the square pore, as we have captured the experimental data with our predictions and computed sensible probability distributions for the bridging time. The model predictions applied to the wave pore are reasonable, giving evidence of the similarities between the cell migration and cell proliferation mechanisms between the two geometries, although this is di fficult to assess with few data points. Section 3.2.2 discusses the analogous results where we instead consider the model predictions from the experimental data on the wave. 3.2.2. Wave-like pore D𝜆 [𝜇m2/day2]0 200400600𝜆 [1/day] 0.02.55.0(a): Likelihood CRs D𝜆 [𝜇m2/day2]0 200400600ℓp(D𝜆) −2−10(b): Profiles for each u0t [day]0 10 20 3 0 40𝜇c(t) 0.00.250.50.751.0(c): Wave coveraget [day]0 10 20 30 40𝜇p(t) 0.00.250.50.751.0(e): Wave perimet ertb [day]12 16 2 0p(tb) 0.00.10.2(g): Wave p (tb) t [day]0 10 20 3 0 40𝜇c(t) 0.00.250.50.751.0(d): Square covera get [day]0 10 20 30 40𝜇p(t) 0.00.250.50.751.0(f): Square perimetertb [day]15 1 923 27p(t b) 0.00.10.2(h): Square p(tb)u0 0.2 0.3 0.4 Figure 8: Likelihood analysis results for the wave pore. In (a), the lines give the boundaries of the 95% credible region for θfor each u0, and the vertical dashed lines show the MLE for Dλ(see Table 1). The profile log-likelihoods for Dλfor each u0are shown in (b), with the threshold c∗≈−1.92 shown with a horizontal red line and the vertical dashed lines show the MLEs for Dλ. In (c)–(f), predictions for µc(t) andµp(t) on each pore geometry are shown, with the blue dots showing the experimental data, the surrounding solid lines giving 95% prediction intervals for each u0, and the dashed lines showing the corresponding estimates at the MLE ˆθ. The estimates for the PDF p(tb) of the bridging time on each pore geometry are shown in (g)–(h). The results in (d), (f), and (h) are predictions on the square geometry using parameters inferred from the wave pore data. We now show the results obtained when we instead estimate Dλfrom the data on the wave pore. The results have greater uncertainty here than in the square case since we have far fewer data points and only one day is covered by the data. The credible regions in Figure 8(a) are much wider and flatter at the bottom than they were in the square case (Figure 5(a)), meaning the estimates for Dλare less precise. The credible region is not well-defined for larger values of Dλ, indicating that we are unable to give any estimate for an upper bound on Dλ. The corresponding profile log-likelihoods in Figure 8(b) for each u0do not intersect the threshold c∗≈−1.92, independently of u0, and so we are unable to give any estimate for the upper limit of the confidence intervals for Dλ, as was already suggested from Figure 8(a). Despite these di fficulties, Figures 8(c)–(f) suggest that we are able to recover values for the experimental data yi,j candyi,j pon each pore geometry, with the predictions from the MLEs going through the experimental data points on each geometry. The probability distribution for the bridging time for the wave pore is shown in Figure 8(g), 15 where we see a mode for tbaround 14 days, with most pores predicted to close between 10 and 18 days, which is consistent with Figure 4. The corresponding probability distribution for the square geometry is shown in Figure 8(h), where we see that most pores are expected to close between 15 and 25 days, which is a shift from Figure 5(g) but is still consistent with the experimental images in Figure 3. It is important to emphasise that the recovery of these results on the square pore with so few data points is remarkable, as it (1) demonstrates the similarity between the migration and proliferation mechanisms on the two geometries, providing further evidence for the first research question, and (2) shows that we do not need such detailed, or even plentiful, data to recover these cellular mechanisms from another geometry. y [𝜇m] 0200400600Day 7Day 1 4Day 2 5Day 28y [𝜇m] 0200400600 x [𝜇m]0200 400600y [𝜇m] 0200400600x [𝜇m]0200400600x [𝜇m]0200 400 600x [𝜇m]0 2 00400600𝜃̂L 𝜃̂ 𝜃̂U Figure 9: Model predictions for the variability in the tissue growth behaviour for the wave geometry for u0∈{0.1,0.2,0.3,0.4}at days 5, 7, 14, 25, and 28. The parameters used are ˆθL=(ˆDL,ˆλL)=(581µm2/day,0.18 day−1),ˆθ=(ˆD,ˆλ)=(58µm2/day,7.29 day−1), and ˆθL=(ˆDU,ˆλU)= (117µm2/day,10 day−1). Predictions of tissue growth for the wave geometry, based on experimental data on the wave geometry, are shown in Figure 9 for u0=0.2. These model simulations are consistent with our experimental observations in Figure 4. The middle row, displaying the most likely outcome, shows that essentially all pores will be closed by 14, which matches Figure 4. The model simulations at day 14 show that while some pores may be completely open at this time, some may be closed or almost closed, and we expect this significant variability as we have so few data points. The predictions we make for the tissue growth on the square geometry using model results on the wave pore are 16 y [𝜇m] 0200400Day 7Day 1 4Day 2 5Day 28y [𝜇m] 0200400 x [𝜇m]0 200 4 00y [𝜇m] 0200400 x [𝜇m]0 200 40 0x [𝜇m]0 2 00 400x [𝜇m]0 200 4 00𝜃̂L 𝜃̂ 𝜃̂UFigure 10: Model predictions for the variability in the tissue growth behaviour for the square geometry using results from the wave geometry for u0=0.2 at days 7, 14, 25, and 28. The parameters used are ˆθL=(ˆDL,ˆλL)=(581µm2/day,0.18 day−1),ˆθ=(ˆD,ˆλ)=(58µm2/day,7.29 day−1), andˆθL=(ˆDU,ˆλU)=(117µm2/day,10 day−1). given in Figure 10. We observe significant variations in the closing time. In particular, while there could be some pores that are only halfway bridged by day 25 or day 28, most are close to closing by day 14 and completely closed by day 25, consistent with Figure 3. The middle row shows that the most likely outcome, according to our model, is that the majority of pores will be halfway bridged at day 14 and closed by day 25, which again matches Figure 3. Overall, these model predictions indicate that, despite only having three data points, we have been able to calibrate our mathematical model su fficiently well so that we capture the original experimental data on the wave from our model predictions and, most importantly, we can recover the experimental data from the square pore from our experimental data on the wave pores. Moreover, the estimated probability distributions for the bridging times on each geometry are a good match to the experimental images in Figures 3–4, as are the predictions of the tissue growth from the model simulations. Thus, not only have we demonstrated the practicality and utility of our method for obtaining these probability distributions, we have provided much stronger evidence than in Section 3.2.1 that the cellular mechanisms driving tissue growth are similar between the two geometries. 17 3.2.3. Tabulated comparison between the two geometries Table 1 compares values of Dλandtb, foru0=0.2, for the two geometries; the choice of u0=0.2 is not significant as the model results are relatively insensitive to this choice. We see that, while we cannot estimate the upper limit of the confidence interval for Dλusing the wave geometry, the lower limits are similar between the two geometries, as are the MLEs. Note that while the MLEs di ffer by a factor of three, this is an insignificant amount when we note that estimates for Dcould vary by many orders of magnitude [15]. These values for Dλprovide strong evidence that the cellular mechanisms driving tissue growth on the two geometries are the same. The estimates of the bridging times for each geometry are not too dissimilar when using either the same geometry or predicting from the other geometry. Square Wave MLE 95% CI MLE 95% CI Dλ[µm2/day2] 152 (90 ,299) 423 (107 ,−) tb(square) [day] 27 (24 ,31) 19 (15 ,22) tb(wave) [day] 21 (19 ,24) 15 (13 ,18) Table 1: Estimates for Dλandtbon each geometry for u0=0.2. The 95% CIcolumn gives the 95% confidence interval for the respective quantities, and the MLE column shows the corresponding MLEs. The second row gives predictions of the bridging time on the square geometry, while the third row is for the bridging time on the wave geometry. Overall, these model results support the hypothesis that the cellular mechanisms driving the tissue growth on each geometry are similar. Moreover, the ability to calculate reasonable estimates and probability distributions for the bridging time provides evidence that the results obtained from one geometry can be used to make predictions about tissue growth on another geometry using the available data. 3.3. Prediction of tissue growth on a hypothetical geometry We now take the results on the square geometry and use them to make predictions on a new geometry that is yet to be experimentally tested. The purpose of this exercise is to demonstrate how our mathematical modelling tools could be used for making predictions on a new geometry from experimental results on a simpler geometry, such as about bridging times, without having to conduct any (potentially expensive and time-consuming) experiments. The geometry we consider is a cross-shaped pore. Using the same values for Dλas were used in making the predictions in Figures 6 and 7, we produce the model predictions for the variability in the tissue growth in this new geometry in Figure 11. We see a similar variance in the previous predictions, namely most pores are closed by day 25 but at ˆθLthere are still some pores that remain half closed. The void maintains the symmetry of the geometry, forming a diamond shape during the early part of the growth process. These results are also largely independent of u0, as we find when plotting these predictions for other u0values (not shown). We similarly show predictions for the yi,j candyi,j p in Figure 12(a)–(b) and the hypothetical bridging time distribution p(tb) in Figure 12(c), all for u0=0.2 An important feature of working with these predictions is that, once the likelihood results have been obtained, producing these predictions in Figures 11–12, or for any new geometry, is not a significantly time consuming task. The snapshots in Figure 11 take only around a minute to compute and visualise, and the model predictions in Figure 18 12 may take around 10 minutes to an hour, depending on the number of samples requested for the prediction intervals. Thus, this type of exploratory analysis of a new geometry can be e fficiently performed in a reasonable time. y [𝜇m] 0200400Day 7Day 1 4Day 2 5Day 28y [𝜇m] 0200400 x [𝜇m]0 200 400y [𝜇m] 0200400x [𝜇m]0 200 400x [𝜇m]0 2 00 400x [𝜇m]0 200 400𝜃̂L 𝜃̂ 𝜃̂U Figure 11: Model predictions for the variability in the tissue growth behaviour for the hypothetical geometry using results from the square geometry foru0=0.2 at days 7, 14, 25, and 28. The parameters used are ˆθL=(ˆDL,ˆλL)=(73µm2/day,1.2 day−1),ˆθ=(ˆD,ˆλ)=(32µm2/day,4.9 day−1), andˆθL=(ˆDU,ˆλU)=(60µm2/day,5 day−1). t [day]010203040𝜇c(t) 0.00.250.50.751.0(a): Void coveraget [day]010203040𝜇p(t) 0.00.250.50.751.0(b): V oid perimetertb [day]11 1417 202326p(t b) 0.00.1(c): p(tb) Figure 12: Predictions for the summary statistics and bridging distributions on the hypothetical geometry using results from the square geometry in Figure 5 at u0=0.2. In (a)–(b), the dashed lines show the estimated curve corresponding to the maximum likelihood estimate for ( Dλ,λ) from the square pore, and the curves surrounding it define a 95% uncertainty band for the curve. 19 4. Discussion New experimental images and modelling predictions in Section 3 provide answers to the three research questions posed in Section 1. The first question, namely whether the cellular mechanisms driving tissue growth are independent of pore shape, appears to be true for the geometries considered. The second question asks whether we can make predictions of bridging times, with uncertainty, on a geometry from a separate geometry. We have found that we can produce reliable predictions with uncertainty between separate geometries, both in the form of probability distributions and prediction intervals. The third question concerns the type and quantity of data required for making predictions between geometries. We found that the data, and the amount of data, we use for summarising the images is su fficient for making predictions of tissue growth with uncertainty, namely information about the tissue void — even with only three data points on the wave geometry. Interestingly, the answers to these questions require only very simple measurements of the experiments, rather than performing cell counting [49]. This observation agrees with previous work that has compared methods using leading edge detection and cell counting, demonstrating that tracking the leading edge is su fficient for estimating the cell migration and cell proliferation rates [49]. We show in Section S6 of the Supplementary Material that if we considered only void area for the analysis of the square geometry then we would obtain the same conclusions, but both area and perimeter are necessary for the wave geometry to answer the research questions. The fact that both area and perimeter are required in general can be expected since area and perimeter together give a detailed description of a shape, but not separately. These answers have important implications. Firstly, we have demonstrated the ability to extrapolate from exper- imental results on one geometry to another geometry, making predictions with uncertainty. This facilitates fast and inexpensive pilot studies to be performed for new pore geometries without conducting the experiments or even fabri- cating the sca ffolds, as with our exploratory analysis in Section 3.3. Numerical simulations could be performed in a few minutes of computation on a standard desktop computer, while conducting the necessary experiments will require more than one month for tissue growth and a considerable amount of e ffort and expense to fabricate the sca ffolds with melt electrowriting. We do not mean to imply that these predictions can replace experimental verification, instead we view this suite of predictive tools as complementary screening tools that can be used to plan and interpret experi- ments e fficiently. Secondly, the novel method we present for obtaining probability distributions for the bridging time provides a useful tool for meeting certain needs by helping us to understand the amount of time required for a tissue to form and bridge the pore, and also for understanding how long tissue growth needs to be incubated for in tissue engineering constructs. Together with the type of exploratory analysis demonstrated in Section 3.3, these probability distributions can help facilitate the construction of a geometry that is likely to bridge within some time window for certain clinical needs, and for determining how long an experiment should be run for. The mathematical model we use in this study is relative simple as it involves just three parameters: D,λ, and u0. A model that better incorporates other e ffects such as cell adhesion or the di fferent phases of tissue growth, in particular the initial phase where cells move o ffthe sca ffold or the later phase where the pore is closing and cells overlap [28], 20 could be of interest to provide more biological insight. A key limitation of working with a more detailed mathematical model, however, would be the need to collect significantly more data so that the necessary additional parameters can be properly estimated [40, 45]. 5. Conclusion and future work In this study, we use a reaction-di ffusion model together with a likelihood-based uncertainty quantification frame- work to study how pore geometry a ffects tissue growth, particularly in how we can make inference about tissue growth on complicated pore geometries using data from tissue growth on simpler square geometries, providing new tools for studying tissue growth with uncertainty and providing probability distributions for bridging times. We use data from pore bridging experiments to perform this analysis, considering a square geometry and a wave-like geometry. Our combined experimental and mathematical modelling results suggest that the cell migration and cell prolifera- tion mechanisms driving tissue growth appear to be independent of the pore geometry, giving evidence that observed curvature e ffects are due to space availability rather than cellular mechanisms. We can make predictions of the bridg- ing time on a new geometry in the form of a probability density function, a powerful tool for understanding both quantitatively and qualitatively what may happen in a pore bridging experiment on a new geometry, including the estimation of probabilities of bridging times over a given time interval. There are several avenues for future work based on our findings in this study. First, our computational tools can be applied to new pore bridging experiments involving di fferent geometries or di fferent cell lines since our methods are independent of these two features. Secondly, it would be of interest to collect more data across more time points to explore the extent to which additional parameters, such as Dandλseparately, can be estimated [15]. If future works consider more than one variable, we note that it would not be feasible to work with the plots of the log- likelihood function as we have done, and instead the profile log-likelihood would be required to obtain uncertainty quantification. Thirdly, the ability to make predictions on new geometries can facilitate a systematic study of how bridging times depend on curvature, such as by defining a geometry that depends directly on a specified curvature and comparing the probability distributions over many curvatures. This would provide a plausible set of outcomes to be analysed prior to running full-scale experiments exploring these features. The predictions made on geometries from data on a separate geometry can be useful for facilitating a pilot study for pore bridging experiments on the geometry, such as the geometry demonstrated in Figure 11 and Figure 12. It would be of interest to see how well these predictions can help with preparing and investing into future experiments, for example in estimating what time scales an experiment may need to be run for by assessing the uncertainty in the bridging times. All code and data to reproduce this work are available on GitHub at https://github.com/ DanielVandH/PoreBridging.jl in the J ulia language [50]. Acknowledgements We thank two anonymous referees for their helpful comments. 21