In-silico and in-vitro morphometric analysis of intestinal organoids

Organoids offer a powerful model to study cellular self-organisation, the growth of specific tissue morphologies in-vitro, and to assess potential medical therapies. However, the intrinsic mechanisms of these systems are not entirely understood yet, which can result in variability of organoids due to differences in culture conditions and basement membrane extracts used. Improving the standardisation of organoid cultures is essential for their implementation in clinical protocols. Developing tools to assess and predict the behaviour of these systems may produce a more robust and standardised biological model to perform accurate clinical studies. Here, we developed an algorithm to automate crypt-like structure counting on intestinal organoids in both in-vitro and in-silico images. In addition, we modified an existing two-dimensional agent-based mathematical model of intestinal organoids to better describe the system physiology, and evaluated its ability to replicate budding structures compared to new experimental data we generated. The crypt-counting algorithm proved useful in approximating the average number of budding structures found in our in-vitro intestinal organoid culture images on days 3 and 7 after seeding. Our changes to the in-silico model maintain the potential to produce simulations that replicate the number of budding structures found on days 5 and 7 of in-vitro data. The present study aims to aid in quantifying key morphological structures and provide a method to compare both in-vitro and in-silico experiments. Our results could be extended later to 3D in-silico models.


Introduction
Organoids are three-dimensional cell cultures capable of self-organisation which can replicate organ-specific shapes, such as crypt-like invaginations in the case of intestinal organoids.Organoids present a powerful tool for studying morphological phenotypes in diverse tissues such as the brain [1][2][3], intestine [4][5][6], and kidney [7,8], among others.However, the mechanisms underlying the generation of heterogeneous morphologies in organoids are not entirely understood, and thus, it is still a challenge to standardise their cultures for experimental manipulation.The development of mathematical models capable of simulating the spatiotemporal dynamics of these complex systems may help to clarify the mechanisms that give rise to the tissue-specific organoid morphology.Mathematical models could also aid in standardising culture protocols, increase experimental efficiency, and lead to improved clinical treatments and regenerative therapies [9,10].Nonetheless, mathematical models require quantitative data to be calibrated and tested.
Due to their self-organisational and heterogeneous properties, the morphological measurement and analysis of organoids can be demanding, particularly if no staining or fluorescence is added to the culture and only bright-field images are obtained.The complexity of the 3D environment within an organoid provides a more faithful recapitulation of in-vivo conditions than other in-vitro approaches.However, the 3D structure can make image analysis of organoids challenging.
To date, organoid morphology has been quantified by segmenting in-vitro images and extracting basic morphometric measurements.Segmentation is a fundamental step in the extraction and analysis of morphological information.Several segmentation algorithms can be implemented on microscopy images, depending on the image quality and other properties.Among those available are region growing, seeded watershed, k-means clustering, and active contours [11].Several segmentation methods have been developed for labelled and label-free 2D and 3D cultures [6,[12][13][14][15].Once an image is segmented, typical morphology metrics such as diameter, perimeter, length, area, circularity, curvature and sphericity can be easily calculated [16].However, there is still a need for algorithms that are capable of aiding in the automatic detection and quantification of prominent tissue-specific features, such as crypt-like domains in intestinal organoids.
The intestinal organoid culture, described by Sato et al. [17], was the first three-dimensional cell culture model capable of self-organising and producing intestinal crypt-and villus-like compartments.The crypt-like units contain the stem-cell niche, while most differentiated cells reside in the villus-like sections.The crypt-villus structures are the most critical characteristics of the intestinal epithelial morphogenesis [18].The presence and frequency of crypt budding is a crucial parameter in the evaluation of organoid-forming efficiency [19].Similarly, this parameter can be used to determine the maturity of an organoid during experimentation [20].Usually, for these cultures, the experimentalist manually quantifies the number of crypts per organoid by observing a small set of organoids or measuring circularity as a proxy to infer crypt growth.
In recent years, there has been a great development of mathematical and computational models of organoid growth [21].Several models have been developed to understand relations between signalling pathways, cell differentiation regulators, and their effect in the growth and shapes of intestinal organoids using different simulation frameworks [22][23][24].However, comparatively few models have focused on possible biomechanical cell interactions involved in the creation of crypt-like structures.Langlands et.al. [25] proposed a model in which crypt fission results from the interaction of two different cell populations with different stiffnesses, which the authors related to stem cells and Paneth cells.This model was developed using the Chaste framework [26,27], which is a popular software library that can be employed to simulate multicellular tissue populations [28,29].Model simulations were compared to in-vitro experiments by measuring the organoid's circularity.Almet et.al. [30] further explored the biomechanical properties involved in crypt fission in this model by performing a parameter sweep over the hard cells' stiffness ratio and their target proportion.They found that these two parameters affect the deformation of the epithelial monolayer and the generation of crypt-like structures.However, in both studies, the model's assumptions for cell proliferation were not based on a specific cell type characterisation.
Similarly to the previous study, they employed circularity to determine crypt fission occurrence in their simulations.
Nevertheless, relying on the circularity measure to indirectly quantify the number of crypts does not accurately quantify organoid growth once they lose their initial spheroid shape.To address this issue, we present an algorithm that uses both segmented in-vitro and in-silico organoid images and retrieves an approximated number of crypts per organoid.Additionally, we propose modifications to the model proposed by Langlands et al. [25] to improve its biological basis and we compare both original and modified model simulations to in-vitro organoid images.Our results should provide insights into the multiscale processes orchestrating intestinal organoids growth.

Segmentation
Brightfield images were acquired on days 3, 5 and 7 of intestinal organoid culture, maintaining the samples at 37˚C, using a Leica DMI6000 inverted microscope with 5x and 10x magnification lenses from the Wolfson Bioimaging Facility, University of Bristol.Images were acquired as 30 z-stacks to cover the depth (z-axis) of the Matrigel domes.First, all images obtained from the microscope were pre-processed by stacking the focal planes in the z-plane to get a single-focus image.This technique combines multiple focal planes into a single image, ensuring a clear focus.We utilised the open-source 'fstack' function [32] to perform this process.The function employs a noise-robust selective all-in-focus algorithm to identify the focused sections from the image groups captured across a focal range.These focused sections are then compiled to generate a high-quality 2D all-in-focus image.Second, two segmentation methods were compared to assess the effect of segmentation quality on the proposed automatic crypt quantification: (i) manual free-hand segmentation and (ii) an open-source segmentation software specially developed for analysing organoid cultures.
Manual segmentation.The manual segmentation was performed in MATLAB using the in-built function imageSegmenter and the assisted freehand to select the regions of interest (ROI).Viable organoids were selected through visual inspection by the user.
Open source-segmentation software for organoids.OrganoSeg is a software developed using the Image Processing Toolbox from MATLAB [33].To use this software, the images (grouped by date) were used as input and the values of intensity threshold (0.1-0.7), max window size (100-500 pixels) and size-exclusion threshold (500-2500 pixels) were selected for each image separately.
The OrganoSeg software presents a few advantages compared to the previous method (i) as it was developed especially for organoid segmentation.One of these advantages is the option to detect all elements that comply with the selected parameters.The user is then given the option to remove detected objects that the user can visually classify as wrong.Another useful option that the interface includes is the ability to select an out-of-focus correction, which can be especially useful for focal plane stacked images.However, using this method makes it difficult to segment organoids that overlap or touch each other.

PLOS COMPUTATIONAL BIOLOGY
All output masks from both methods were saved for later analysis.The output masks contained those segmented organoids found to be viable.The user classified whether each organoid was viable by visually considering if the membrane on the organoid appeared intact and there were no breakage and leak of cells into the Matrigel.

In-silico model
Langlands et al. model [25].The first computational model that explored the effect of heterogeneous biomechanical cell properties on the development of crypts was performed by Pin et al. [34].Later, a study by Langlands et al. [25] investigated this effect on a two-dimensional computational model of intestinal organoids.This model was then extended by Almet et al. [30] to explore further biomechanical properties involved in crypt fission.The general premise of these 2D cell-based models was to test if crypt fission could be initiated solely by different stiffness of two cell types, one soft and another hard, rather than by a signalling cue.
The model in [25] represents a confluent epithelial monolayer (i.e. a cross-section of a 3D intestinal organoid) and was developed using the agent-based computational framework Chaste [26,35].The monolayer is defined by a set of epithelial cells delineated by a honeycomb mesh created using Voronoi Tessellation and Delaunay triangulation.In the system, cells become neighbours and connected if they share an edge in the Delaunay triangulation.A cell can lose contact with another one if a set cut-off length distance is exceeded.The epithelial layer is formed by two cell types, which differ only in their relative stiffness.During the simulation, cells are exposed to interactive forces between their neighbours and the simulated Matrigel through a linear spring force.For more information about parameter values please see Table 1.
The proliferation of all epithelial cells in this model is defined by a stochastic cell cycle duration sampled from a uniform distribution of U (12,14) hours.When a cell divides, two daughter cells connected by a spring laying in a randomly chosen direction instantly replace the original cell.The daughter cells are initially positioned in opposing orientations, 0.05 cell The majority of parameter values are maintained between the original Langlands et al. [25] and our modified unless stated otherwise.Distances are measured by cell diameters (CD), and time is scaled in hours. .Thus, both cell types will always be present in the system according to the set proportion.Additionally, all epithelial cells divide ad infinitum.
Cell death is defined by anoikis: any epithelial cell that loses contact with its neighbours and enters the lumen or the Matrigel is removed from the simulation.Otherwise, if the epithelial cell is not separated from the layer, it can keep dividing according to the set proportion.This model has proved capable of mimicking crypt formation, showcasing the effect of biomechanical properties, like stiffness and the cellular target proportions of soft and hard cells required to promote crypt fission using only two cell populations [25,30].However, there is still room for improvement, as its assumptions for cell proliferation are not based on a specific cell type characterisation.The first study that describes the model [25] refers to the model by Pin et al. [34], in which the Young modulus (stiffness) of Paneth cells was approximated in comparison to that of stem cells.Their results showed that Paneth cells are approximately four times stiffer than stem cells.Thus, it could be assumed that stem and Paneth cells are soft and hard cells in the model, respectively.However, the model does not take into account the presence of other intestinal cells such as enteroendocrine (EEC), enterocyte (EnC), transit-amplifying cells (TA), goblet (GC) and Tuft (TC) cells.Additionally, it is known that Paneth cells are terminally differentiated cells and not capable of producing stem cells nor other Paneth cells [39,40].Therefore, with the current 2D model configuration it would be challenging to produce crypt fission if Paneth cells were the most stiff cells in the system, as other cells would not be capable to push or move them.
Therefore, we propose a few modifications that aim to enhance the model's representation of biological basis such as the initial cell population proportions and the model's cell proliferation.These modifications are further explained in the Results section below.

Crypt counting method
The organoids' masks must be processed to obtain possible crypt sections.S1 and S2 Figs summarise the crypt counting method processes applied to in-vitro and in-silico organoid images respectively.One important difference between the in-vitro and in-silico data was the method of extraction used.While in-vitro organoid images are segmented to obtain the organoid ' In principle, one could apply either boundary extraction mechanism to in-silico data.The cell-centroid method was less computationally intensive and minimised user input, given that it automatically returned information on the organoid boundary.However, one can use any other method to extract the boundary and optimise the parameters accordingly, following the method described in the training section below.
Each raw boundary is estimated using a Fourier approximation [42] to soften the mask's rough edges acquired during the segmentation image processing (see S1(d) and S2(d) Figs).It is expected that the required number of harmonics used to approximate the boundaries would change between experimental days as organoids grow and have more crypts with less definition reflected in the mask.In the images obtained from three-day-old organoids, it can be observed that most organoids have between 1 and 2 crypts, and the arc length of the crypts is relatively large compared to the main body of the organoid.Therefore, these images do not require high precision, and a lower number of harmonics suffice, as it produces a softer boundary approximation.However, as organoids grow and produce more crypts, these cannot be defined as clearly in the obtained masks, as can be observed in  data for days 5 and 7.In these cases, small changes in curvatures can indicate the presence of a crypt, and the Fourier approximation is required to be more precise.Thus, images obtained on days 5 and 7 of the experiment would need a higher number of harmonics.Therefore, such number was included as a main parameters of the algorithm used to detect crypts and required optimisation.These hypotheses are confirmed in the values obtained by the optimisation shown in Table 2.
\In the present study, we used the approximated boundary and calculated its curvature and normal as functions of the boundary location.The curvature and normal values were then used to modify the boundary and obtain all convex and concave regions (see S1(e) and S2(e) Figs).A middle point was calculated for all concave regions (inward sections), and these points were used to select regions and decide if a crypt was present on that region (see S1(f) and S2(f) Figs).
Parameterisation.We defined a set of main parameters that helped defining a crypt-like section automatically.These parameters were the number of harmonics required for the boundary approximation, the potential budding section area and arc length.However, for the last two parameters, we decided to consider normalised values according to the total organoid value, as these will vary according to the size of organoid image: Budding section area Total organoid area ; ð1Þ ArcLength norm ¼ Budding section arcLength Total organoid arcLength : ð2Þ Our approach therefore accounts for the expected budding size range according to the experimental day.Each day should have its optimal value of the minimum and maximum normalised area (Area norm ), allowing us to discriminate noise and villi regions from crypt regions.Additionally, a threshold for minimum normalised arc length (ArcLength norm ) was included to avoid areas with sharp corners that are not present in crypts.
Training.Our experimental data consists of a total of 69 images: 30 of day 3, 19 of day 5, and 20 of day 7.The selection of best-fitted parameters was calculated by training the code using a selection of 5 images of each experimental day, which contained a different number of organoids (i.e.day 3 data contained 18 organoids) in the case of in-vitro data.In the case of insilico data, a group of 10 organoids per simulated day was selected for training.All the selected organoids' crypts were manually counted as a control during the training phase.
An objective function was generated for the training and parameter optimisation.This function included the code-calculated crypts, the hand-counted crypt values, and the calculation of percentage error comparing these two values: This objective function was used to find the global minima using the simulated annealing algorithm within MATLAB's Global Optimisation Toolbox.This method was selected due to its capacity to explore globally for parameter values and prevent the system from being trapped in local minima in early iterations.
The optimised parameters obtained from the training are displayed in Table 2.In general, the training data suggest that the code can have around *80% of accuracy, so it can be used to approximate the mean crypt count at different time points of an intestinal organoid experiment.
The optimised parameters were applied to calculate the number of crypts on the manually segmented organoids (n = 132 organoids per day) and the organoids obtained with the OrganoSeg algorithm (n = 124 for day 3, n = 122 for day 5, n = 109 for day 7).Additionally, as a control, we counted the crypts from the manually segmented organoids.For the in-silico data, on both the 'original model' and the 'proposed model', we had 50 simulated organoids per experimental day.

Circularity
Previous studies have used circularity as an indirect method to measure the growth of intestinal organoids [25,30,[43][44][45].Circularity has been regarded as an effective method due to the initial spheroid shape of organoids and their subsequent production of crypts at later stages produces a loss of "roundness" or circularity.Therefore, it can be expected to obtain an almost perfect circle when the organoid has not yet produced crypts (circularity ' 1) and to lose this roundness once crypts appear (circularity <1).Here, we will compare the results obtained from the presented crypt count method and the circularity obtained from the same images.Circularity can be calculated as: All p-values shown in Figs 3, 6, S4, S5 and S6 are between each indicated set of bars, and were calculated using a two-sample t-test.

Results
Fig 3 exhibits the results of the mean number of crypts obtained from manual user count ('usercounted') and the code-counted using the previously presented segmentation methods, namely, OrganoSeg and Manual segmentation.To avoid confusion we will name the results obtained

PLOS COMPUTATIONAL BIOLOGY
from applying the crypt-count algorithm to masks obtained from OrganoSeg and Manual segmentation 'OS-Code-count' and 'MS-Code-count' respectively.Here, it can be observed that the code-counted results obtained from 'OS-Code-count' and 'MS-Code-count' are not significantly different at day three compared to 'user-counted' data.The overall trend among the three data sets shows similar organoid growth.However, in data sets where the code-counting crypt algorithm was applied, we can observe a more significant increase in the mean number of crypts between days 3 and 5; and a minor growth change between days 5 and 7.
Results obtained from 'MS-Code-count' show no significant difference on days 3 and 7 compared to 'user-counted' data.Nonetheless, results from day 5 show a substantial difference between these two sets.S3 These results suggest that segmentation accuracy is essential to improve the performance of the code-counting crypt program.In the case of 'MS-Code-count', possible factors that limit the accuracy at day 5 include the variability in crypt sizes or the inability of the optimisation function to obtain better parameters at this stage.As mentioned in the methods section, the code is limited to a specific set and parameter range that may not be as flexible to cover the variability in organoid morphology present at this point of the experiment.
Next, we applied the crypt-counting algorithm to the in-silico organoids obtained using the Langlands et al. model [25].

PLOS COMPUTATIONAL BIOLOGY
replicating organoid growth in terms of crypt count at day 7.However, there was room for improvement to replicate the growth rate observed on days 3 and 5. Therefore, we decided to perform some modifications to the model to make it more biologically realistic and observe if these modifications were enough to improve the growth rate of the number of crypts in the insilico organoids.

Proposed modifications to the in-silico model
Initial conditions.As previously mentioned, the in-silico organoids are compared to the control in-vitro experiment using crypts extracted from mouse intestine.Several studies mention that the intestinal organoids contain similar cell proportions that those found in-vivo [4,17,[46][47][48][49].A study by Haber et al. [50] reported the profiling of 53,193 individual intestinal cells in the small intestine of 7-to 10-week-old mice.The authors classified and quantified the presence of several cell types as enteroendocrine (EEC), enterocyte (EnC), enterocyte progenitor (EnC-P), transit-amplifying cells (TA) (including on S and M cell-cycle phases), Paneth (PC), stem (SC), goblet (GC) and Tuft (TC) cells.Here, we extracted the cell population numbers of only the stem, TA and Paneth cells found in this study and inferred the initial proportion for the model.We assume that only those cell types are present in the system at the start of the simulation, as the experiments initiate from crypt fractions, and those cell populations mainly reside there.Thus, the initial cell population proportions used on all simulations of the 'proposed model' were: 42% stem cells, 49% TA cells and 9% Paneth cells.On the other hand, the initial proportions set for the simulations performed with the initial model were directly related to the target proportion set for cell production too, which in this case was *20% stem cells and *80% Paneth cells.
We decided to address the limitations we observed in the initial model and performed a few modifications to observe their effects on the resulting simulations.First, the initial number of cell-nodes in the system was increased from 20x20 to 36x40, to allow greater area for the organoid's growth during more extensive simulations.Additionally, the initial organoid centre was moved to fit the centre of the new mesh.
Cell production model modifications.In the previous model, a hard cell could divide and produce a soft cell, and vice-versa.However, if we associate stiffness with a specific cell type, their division would be restricted to their lineage.As previously mentioned, stem (SCs) and Paneth cells (PCs) have been previously represented as soft and hard (i.e., stiff) cells, respectively.Nevertheless, Paneth cells are terminally differentiated and do not divide [51].Suppose the cell production relies only on stem cells in the simulation with time.In that case, the system will be saturated by the stiffer cells that cannot divide, as soft cells will be expelled from the epithelial monolayer.
Therefore, we decided to introduce two new cell types.The first was defined as transitamplifying cells (TAs) to represent the possibility of having a highly replicative cell type as a hard cell.A second cell type was defined as a general differentiated enteroid cell (ECs) that represents all other differentiated cells that TAs can produce.We performed a few tests (see S4 and S5 Figs) and found that the simulations that could replicate in-vitro crypt-count numbers more closely were obtained using a 4.5 greater stiffness ratio than SCs in all non-stem cells (i.e., PC, TA, and ECs).
To keep our model simple, we did not include signalling processes that regulate the differentiation of proliferative cell populations.Instead, our model keeps relying on a stochastic process and suggested cell proportions to generate the different cell populations.Stem cells in the intestinal crypt have been extensively studied and are known to retain their stemness [52,53].According to this existing understanding, the simulated cell cycle was modified to allow stem cells to symmetrically divide and generate other stem cells, which simulates the maintenance of multipotency; or to asymmetrically divide and produce a TA or Paneth cell daughter (see Fig 2b).The stem cell division will be stochastic, depending on a specific target proportion set.In this study, the probability proportions for the production of each cell type derived from stem cells were assumed to be: • * 89% other stem cell; • * 9% a Paneth cell; • * 2% a TA cell daughter.
Using these proportions, we could regulate the cell production in the model.The mentioned proportions were selected as an educated guess, as PCs consist of around *9% of the cells in the intestine [50] and the model assumes that SCs are their only source.Furthermore, the probability of SC symmetric division (SC prop) was set to (* 89%) based on a study performed by Itzkovitz et al. [54], in which they observed a high symmetrical division of Lgr5+ cells on small crypts.On the other hand, the proportion of TAs obtained from SCs (TA prop) was set to be low (* 2%) to avoid overpopulation of this cell type, as TAs already in the system will produce more TAs by symmetrical division.It is widely known that TAs go through several rounds of cell division before reaching maturation and differentiation [55,56].Studies suggest that TA cells divide between two to five times before committing to a differentiated cell type [57][58][59][60].In our model, TA cells can divide symmetrically.And, to account for the fact that TA cells have a higher chance of differentiating after several cycles, we increased the proportion of differentiated cells produced by TA cells after day 5.This change enables several rounds of cell division to occur before differentiation, without eliminating the possibility of TA cells differentiating before day 5, as each cell starts at a different point in its cell cycle during simulations.
Thus, we assumed that during the first five simulated days of organoid growth, there is a higher probability of TAs producing more TA daughters.Then, the probability of TAs producing ECs is increased after this period.ECs, as well as PCs, are considered final differentiated cells and do not divide (see Fig 2b).
The probability proportions for TAs cell division are assumed to be: • Before day 5 • * 90% other TA daughter • * 10% ECs daughter • After day 5 • * 70% other TA • * 30% ECs daughter A sweep of TAs and ECs probabilities was performed and the values presented in this study resulted in the most similar crypt count results compared to our in-vitro experiments.Fig 5 displays images of the simulated organoids using either the original (i.e.[25,30]) and our 'proposed model'.
Results obtained from the mean number of crypts from the 'user-counted' show an exponential increment as the experiment progresses.This effect is displayed in the in-silico results obtained using the 'proposed model'.Fig 6 presents a comparison of the 'user-counted', the 'MS-Code-count' and the code-counted crypts obtained from the 'proposed model'.The results obtained from simulations on day 3 present a higher average number of crypts compared to the number of crypts found in the experiment, similar to what was observed previously in the 'original model'.Nonetheless, for day 5, we observed that the 'proposed model' was able to replicate the mean number of crypts as the in-vitro data quantified by the user.Similarly to what was observed in the initial model, results obtained from day 7 present a similar average number of crypts per organoid across in-vitro sets.
Finally, we were interested in observing how our crypt-count results compare to traditionally used methods.Table 3 and S6 Fig. show the results obtained from measuring the circularity of organoids through manual segmentation on experimental data, and on simulated data (from the Original [25,30] or our model).We should expect reduced circularity with an increased number of crypts as time progresses.Instead, the mean circularity obtained from manually segmented images of in-vitro organoids shows an increment, which suggests greater roundness, from day 5 to day 7. Similarly, the simulations obtained from the 'proposed model' do not show a difference in circularity as organoids grow.

Discussion
Cell and tissue imaging and image analysis can be instrumental to measure and quantify complex cell phenotypes in time and space.There has been a rise in the development of automatic and semi-automatic algorithms that specialise in the segmentation of two and three-dimensional cell cultures, including those based on label-free images [13,[61][62][63][64].
In the case of organoid cultures, algorithms have recently been proposed to segment and perform basic morphological image analysis by using traditional segmentation methods [33], or based on organoid staining [65]; other methods leverage machine-learning/deep learning approaches [66][67][68][69].Nevertheless, these software are yet to be improved to detect and measure  Results obtained from in-vitro and in-silico data obtained from both models.The circularity was calculated using Eq (4), and here the manually segmented images were used.The mean number of crypts per organoid was obtained by applying the counting-crypts code to the resulting simulations of both models (original and proposed).
The mean number of crypts calculated using the models is being compared to in-vitro user counted crypts.The standard error of the mean (SEM) was calculated for each set. https://doi.org/10.1371/journal.pcbi.1011386.t003 unique structures for specific organoid tissues.For example, in the case of intestinal organoids, the appearance and number of buds that replicate the shape and function of intestinal crypts.Previous studies have counted crypts or buds on intestinal organoids by staining or labelling stem cells and segmenting regions containing a group of these cells to automatically or manually count them [70][71][72][73][74].In cases where no labelling or staining is implemented during experimentation, we still rely on the user to define and count the number of crypts per organoid.This occurrence can be attributed to the representation of the initial buckle in the tissue, which may resemble crypt-like structures due to the 2-dimensional nature of our model.Thus, we believe there is an opportunity for new software targeted at quantifying morphological structures on label-free bright-field images of organoids.This measure, especially on intestinal organoids, can be essential to define organoid cultures' proper growth and standardisation.
The crypt-counting algorithm developed in this study presents a valuable tool to calculate an average of crypt-like or buds structures present in their in-vitro cultures.The algorithm has proved capable of approximating the average budding structures found in in-vitro intestinal organoid cultures on days three and seven after seeding.We believe that our method is applicable to other types of organoids that exhibit tissue buckling and branching detected in 2D images.This includes organoids from various tissues, such as lung, mammary glands and kidney organoids, for example.Our algorithm can be trained to identify finger-like structures of varying sizes and shapes, making it adaptable to different types of organoids.
Nonetheless, the code shows limitations in calculating the number of crypt-like structures in sample dates in which there is a higher variability of organoid growth in terms of the number of budding sections.As the in-vitro experimental data available for our current study was limited, we were unable to incorporate growth-tracking data for the organoids.However, future studies utilising organoid tracking data have the potential to offer a more accurate and comprehensive morphometric analysis.
Additionally, in the case of the crypt-counting algorithm applied to in-silico images, we can observe a mismatch in the mean number of crypts compared to the in-vitro data.On day three, organoids start breaking symmetry from their original spheroid shape.Thus, it could be challenging to distinguish a crypt forming in the early stages, which can lead to an under-or over-counting of crypts, as the training depends on the user expertise.Also, the versatility of the algorithm to be applied to agent-based simulations of organoids provides the opportunity to improve such models and make predictions of pattern or structure formations before experimentation, which has been tested before for pattern formation in 2D stem-cell cultures [75].Due to the limited amount of data available to us at the moment of this study, we were not able to develop a deep learning algorithm to automate crypt counting.However, we believe that in future work, the code limitations could be improved by implementing deep learning and machine learning algorithms to aid in finding additional or better-fitted parameters to define crypt-like structures [67,76,77].
It is important to notice that circularity was used before to test the accuracy of the 'original model'.However, the results obtained in this study are not the same as those presented on [25].Their circularity values calculated for 100 simulated hours are lower than those obtained in the present study.This may be due to the different methods used to extract and approximate the boundary from the simulations.Here, we performed the same boundary analysis to both in-vitro and in-silico images, in contrast to calculating the simulated epithelial layer circularity directly from the model's cell connections.Additionally, the 'original model' did not present issues of cell competence or the removal of soft cells from the simulation due to their cell production model, in which any cell can produce a daughter of any other cell type according to the target proportion set by the user.
The proposed changes to the in-silico model maintain the potential to produce simulations that replicate the number of budding structures found on in-vitro experimentation.Nevertheless, the model still presents limitations, such as the lack of signalling pathways that regulate cell production [78], and it is limited to the 2D structure of the epithelial layer.The initial cell population proportions used at the beginning of our simulations were inferred from a scRNA-seq survey [50].However, it has been reported that dissociation and cell enrichment protocols used in scRNA-seq can possibly bias cell type numbers, making direct estimation challenging [79].While we recognise this limitation in our assumptions, in future work, we aim to incorporate methods such as pseudotime analysis [80] and RNA velocity [81] to more accurately estimate the transition probabilities between identified cell states and improve the reliability of our model.
The spatial restriction affects the position of new cells added to the system during the simulations, as it restricts the free movement of cells and confines their mobility to interactions with neighbouring cells.Thus, we were not able to recreate an accurate location of the celltype populations found on small-intestine crypts, in which SCs and PCs are positioned at the bottom.Similarly, although the results obtained from our model indicate that differential stiffness in the system enhances its ability to initiate crypt fission, the framework structure gives rise to the expulsion of softer cells more easily than hard cells, which may prompt the model to end with fewer soft cells (i.e.Stem cells) at the end of the simulation.In this study, we deliberately focused on 2D images to achieve our main objective, which was to develop a method that facilitates direct comparisons between our in-silico and in-vitro data.While we acknowledge the limitations of representing a 3D system in 2D, we firmly believe that analysing 2D images can still yield valuable insights into the system's dynamics and reveal previously undefined aspects.Although it does not capture the complete complexity of a 3D system, studying 2D representations can contribute to our understanding of the underlying processes.Moreover, our ongoing research is dedicated to developing a 3D in-silico model and comparing it with 3D in-vitro images of intestinal organoids.We recognise the significance of investigating the system in a more accurate 3D representation and intend to address this in our future work.

Conclusion
The present study is limited to snapshots at three different days, with no tracking of organoids, which makes it difficult to observe their growth directly.Nevertheless, the present crypt-count code and in-silico model provide valuable tools with which to explore the effect of crypt development on intestinal organoids and stepping stones for future approaches to classifying morphological structure in organoids.It is important to note that the results presented in this paper merely highlight the complexity of the crypt fission phenomenon and suggest that additional factors beyond cell stiffness and subpopulation proportions may contribute to its initiation.Further investigations are necessary to gain a comprehensive understanding of the underlying mechanisms and refine our model accordingly.We do foresee an increasing adoption of agent-based models also in engineering biology, with the aim at designing novel 3D cellular structures and tissues, possibly integrating a description of biomechanical interactions with detailed subcellular whole-cell processes, and adding feedback control to steer phenotypes of interest [82][83][84][85][86].
Fig 1 presents sample images of the organoid cultures obtained on experimental days 3, 5, and 7 after seeding.

Fig 1 .
Fig 1. Timeline images of in-vitro intestinal organoid culture.Samples of in-vitro intestinal organoid culture stacked images with segmented organoid boundaries highlighted in green.The images were obtained at days 3 (a), 5 (b) and 7 (c) after seeding, using 10x (a) and 5x (b and c) objective.Scale bar: 200 μm.https://doi.org/10.1371/journal.pcbi.1011386.g001 https://doi.org/10.1371/journal.pcbi.1011386.t001diameters distant from the original parent cell, and the new spring connection increases linearly from 0.1 to 1 during the first hour of cell-cycle.A user-defined division parameter controls the probability of soft cell self-renewal in the model.All epithelial cells generate one daughter cell that is identical to the parent and a new cell with probability p of being a soft cell and probability 1-p of being a hard cell (see Fig 2a)

Fig 2 .
Fig 2. Cell cycle diagrams.(a) Original cell-cycle model used by Langlands et al. and Almet et al. [25, 30], (b) 'proposed model' modification that requires a production probability for each cell type derived from stem cells.Both models maintain the stochasticity of the system, as the selection of cell type is dependent on a random number generated at each time-step.Cell types: stem cell (SC); Paneth cell (PC); transit amplifying cell (TA); general differentiated enteroid cell (EC).https://doi.org/10.1371/journal.pcbi.1011386.g002 s boundary (see S1(b) Fig) and the boundary of each organoid was extracted using bwboundaries function of MATLAB (see S1(c) Fig), in-silico organoids' boundaries are obtained by collecting the cells' centre points (see S2(b) Fig), and applying a genetic algorithm of the travelling salesman problem [41] to obtain the organoid's boundary (see S2(c) Fig).
Fig 1 of our experimental

Fig 6 .
Fig 6.Mean number of crypts: Modified in-silico model.Comparison of mean number of crypts obtained using the counting-crypts code applied to the in-silico organoids simulated with the 'proposed model', and to manually segmented in-vitro data ('MS-Code-count'); 'user-counted' crypt numbers are also shown.P-values from two-tailed unpaired t-test computed over the data sets shown, *p<0.05,**p<0.01,***p<0.001,****p<0.0001.https://doi.org/10.1371/journal.pcbi.1011386.g006 the mask's information.(a) Initial stacked image of a 3-day-old organoid culture; (b) manual segmentation of an organoid; (c) raw boundary extracted from the segmented organoid; (d) boundary approximation obtained using a Fourier approximation; (e) calculated concave sections (in red) on the boundary; and (f) possible crypt sections detected by our algorithm.(TIFF) S2 Fig. Process performed on in-silico images.This image shows the general processing performed for each in-silico organoid's boundary to calculate the number of crypts presented after extracting the mask's information.(a) Initial image obtained from a 3 day old organoid simulation with several nodes representing stem cells (blue), transit amplifying cells (yellow), Paneth cells (green), Matrigel (pink), lumen (dark blue) and the simulation boundary (red); (b) location of epithelial cells (SC (blue), TA (yellow), PC (green)) extracted from the simulation; (c) raw boundary of the simulated organoid; (d) boundary approximation obtained using a Fourier approximation; (e) calculated concave sections (in red) on the boundary; and (f) possible crypt sections detected by our algorithm.(TIFF) S3 Fig. Mean number of crypts found from in-vitro images: 'user-counted' vs 'MS-Codecount'.Comparison of number of crypts found per organoid at each day by the user in comparison to those found in 'MS-Code-count'.(a-c) Histograms of the distribution of crypts found by the user (pink), compared to 'MS-Code-count' (blue).(d-f) Boxplots of the number of crypts found using the previously mentioned methods, in which the boundaries of the box represent the 25 th and 75 th percentiles respectively of the median (red line) number of crypts found per day, the width of the notch represents a 95% confidence interval around the median, the whiskers extend the most extreme data points, and the outliers are represented individually (red plus signs).(TIFF) S4 Fig. Mean number of crypts: Comparison of in-silico model with different stiffness ratios for TA cells.Comparison of mean number of crypts obtained using the counting-crypts code applied to the in-silico organoids (n = 10) simulated with different stiffness ratios for TA cells: 2x, 3x and 4.5x relative to the spring stiffness of the 'Proposed model'.P-values from twotailed unpaired t-test computed over the data sets shown, *p<0.05,**p<0.01,***p<0.001,****p<0.0001.(TIFF) S5 Fig. Mean number of crypts: Comparison of in-silico model with different SC cell cycle lengths.Comparison of the mean number of crypts obtained using the counting-crypts code applied to the in-silico organoids simulated with the 'Proposed model' that simulates SCs cell cycle between 22 and 24 hrs (n = 50), and the model sampling the cell cycle length of SC between 46 and 48 hrs ('SC cycle 46-48hrs', n = 10), 'User-counted' crypt numbers are also shown for reference.P-values from two-tailed unpaired t-test computed over the data sets shown were not significant, which suggest that our current model is not statistically significantly sensitive to changes in the cell cycle length of SCs and that this parameter does not have a substantial impact on the overall results of crypt-like structures formation in the model.(TIFF) S6 Fig. Mean circularity.Comparison of mean circularity calculated at each day for in-vitro (using the manually segmented images) and in-silico organoids.P-values from two-tailed unpaired t-test computed over the data sets shown, *p<0.05,**p<0.01,***p<0.001,

Table 2 . Optimised parameter values. Simulated day Parameters Error Harmonics Min crypt area Max crypt area Min arc length
Results obtained during optimisation of parameter values to identify crypts on both in-vitro and in-silico organoids.https://doi.org/10.1371/journal.pcbi.1011386.t002