The classical automated creation of tissue-specific meshes involves three steps (see Figure 1 and compare (Bennetts et al., 2014)): 1) image segmentation and ellipsoid fitting, 2) defining geometric input files, and 3) generating geometry and mesh.

The geometry and distribution of individual cells can be determined directly from the processed images of specimen-specific tissues. To obtain a geometric description of the cells, we fitted ellipsoids to a contiguous set of voxels associated with each cell using ImageJ v.1.53f51 (Schneider et al., 2012) ¹ . Eventually, we obtained the cells’ location, radii, and orientation in the tissue. Anisotropic voxel sizes in the 3D images exacerbated the separation of adjacent cells in the z-direction. Therefore, we have resampled the 3D images to obtain image stacks with isotropic voxels. In addition, various filters were employed to reduce noise in the image, enhance contrast, and remove unevenly illuminated backgrounds while preserving low-intensity fine detail. Subsequently, for the sake of reproducibility, an auto-threshold technique was utilised to segment images into cell interior and exterior instead of manually tweaking them. Subsequently, gigantic artifactual objects with an area surpassing 300 μm² were removed from the image slice. Miniature noises and artifacts cannot be filtered out based on this criterion because parts of the cells could be excluded from individual image slices. As an alternative, a volume filter from BoneJ2 v.7.0.13 (Domander et al., 2021) was employed to remove various tiny noises with a volume of less than 200 µm³. Afterwards, holes were filled and touching cells were separated by the Distance Transform Watershed algorithm (Legland et al., 2016). Lastly, the cell shape was assumed to be ellipsoidal and the geometric description of the cells was obtained via 3D ellipsoid fitting (Ollion et al., 2013). The whole procedure was wrapped into an ImageJ script to batch process all images.

To define geometric input files, we utilised a format consistent with previous research to store the characteristic dimensions of the ECM and ellipsoids (Bennetts et al., 2014). Due to the fact that each image was captured for only one particular cartilage zone, a geometrical model incorporating multiple cartilage zones can be constructed.

Numerical models such as the finite element method (FEM) require a discretised geometry. For that, we employed the mesh generator Netgen (Schöberl, 1997) to automatically construct the geometry. From the geometrical parameter in the input text file, a solid model is assembled using the Open Cascade Technology (OCCT) ² geometry kernel, provided in the Netgen Python library. Boolean unions of touching objects were implemented to obtain a single cell volume. We meshed the cellular geometries using a pre-determined mesh hypothesis. The meshes were created in such a manner that their surface and volume error were less than 1% with regard to the cell surface area and volume.

In this work, we suggest the following approach to build a numerical model in an automated manner (see Figure 1): 1) 3D cell segmentation using an artificial intelligence approach on the original images, 2) geometry creation directly from the segmented image using the Standard Triangle Language (STL), 3) conversion to a computer-aided design (CAD) model and meshing.

Before finding a suitable approach, we compared various methods, namely, the Trainable Weka Segmentation (TWS) method (Arganda-Carreras et al., 2017), DeepImageJ (Gómez-de Mariscal et al., 2021), Biomedisa (Lösel et al., 2020) and ZeroCostDL4Mic (von Chamier et al., 2021). Unfortunately, we could not achieve satisfactory results using these methods. A detailed description of the procedure that we followed for each approach is reported in the Supplementary Material.

We chose Stardist because it was developed with a focus on blob-shaped cells in a crowded environment (Weigert et al., 2020). StarDist uses a slightly modified variant of 3D U-Net as a neural network backbone to predict the radial distances (star-convex polyhedron representation) and object probabilities (indication of how likely a pixel is part of an object). There are only a few hyperparameters in StarDist to be tuned. We trained the model on three different platforms to assess hardware configuration requirements for deep learning model training. A workstation with Intel Core i7-10700 (64 GB RAM) without GPU support, and the post-processing node of the HAUMEA HPC cluster at the University of Rostock with a GeForce179 RTX 2080 Ti GPU (11 GB VRAM) and an Intel Xeon Gold 6234 CPU (8 cores, 16 threads) with 1.5 TB RAM were utilised for training the deep learning model.

Because the model input is the original image, no pre-processing step is required. We used the fitted-ellipsoid images as ground truth. These images were generated by fitting ellipsoids to the voxels of the cells using the 3D ellipsoid fitting plugin in ImageJ. This process was performed after applying the Distance transform watershed. The fitted ellipsoid images are represented as 16-bit images, where each ellipsoid is assigned a unique voxel intensity value. These intensity values serve to distinguish individual cells within the image. An example of a slice in a training pair can be seen in Supplementary Figure S10. To increase the diversity of the training data and improve the models’ performance, we applied a basic data augmentation comprising a lateral flip, rotation, and random intensity scaling. We randomly selected 23 images for training/validation and four images for testing the deep learning model. The test images were labelled according to their source, sample number and cartilage zone, such as animal 1 sample 1 top zone (A1S1), animal 2 sample 1 top zone (A2S1), animal 4 sample 1 deep zone (A4S1) and animal 5 sample 2 middle zone (A5S2). The testing data set includes one good-quality image (A1S1), two noisy images (A2S1 and A5S2) and one very noisy image (A4S1). As a precaution against memory issues, the training/validation dataset has been cropped into 92 equal XY-size images of 512 × 512 pixels (78 for training and 14 for validating). In total, we have 6 GB of data.

To assess the model performance, we used the accuracy estimate introduced by Weigert et al. (2020) as a c c u r a c y τ = T P T P + F N + F P (1) where τ is the overlap threshold, true positives TP are the number of correct predictions, unmatched predicted instances are FP false positives, and unmatched groundtruth instances are identified as FN false negatives. A predicted object I _pred and ground-truth object I _gt are considered as one correct prediction if they have an intersection over union (IoU) that fulfils I o U = I pred ∩ I gt I pred ∪ I gt ≥ τ . (2)

We also report precision (see Eq. 3), recall (see Eq. 4) and F1 − score, also known as the Dice score (see Eq. 5). They are commonly employed in practice to access the model performance (Laine et al., 2021). Precision measures the proportion of predicted voxels in the segmentation results that match the ground truth voxels. Over-segmentation results in lower precision scores. On the other hand, recall represents the proportion of predicted voxels in the ground truth that were correctly identified. Low recall scores are caused by under-segmentation. The F1—score, which combines precision and recall, provides an overall assessment of segmentation performance. An increase in F1, indicates better segmentation performance, with a value closer to 1 indicating more accurate and precise segmentation. p r e c i s i o n = T P T P + F P (3) r e c a l l = T P T P + F N (4) F 1 = 2 × p r e c i s i o n × r e c a l l p r e c i s i o n + r e c a l l (5)

Apart from building the model from scratch, we also employed transfer learning techniques with the StarDist model. Our approach involved retraining the pre-trained StarDist model from Rasse et al. (2020), ³ with our specific data. The pre-trained model had been previously trained on an Arabidopsis thaliana lateral root nuclei dataset (Wolny et al., 2020). Furthermore, for comparison, we considered ilastik (Berg et al., 2019) and Cellpose3D (Eschweiler et al., 2022). A comprehensive description of the methodology employed in implementing ilastik and Cellpose3D can be found in the Supplementary Material.

For the geometry creation, we generated the STL mesh of the cells directly from the deep learning predicted images rather than fitting the ellipsoid to derive the characteristic dimensions of the cells. scikit-image ⁴ was leveraged to accomplish this task using the marching cubes algorithm (Lewiner et al., 2003). To improve the STL mesh quality, an original surface representation was remeshed and smoothened employing isotropic explicit remeshing (Peercy, 1993) and Taubin smoothing (Taubin, 1995) available in PyMeshLab (Muntoni and Cignoni, 2021). A simplification of the resulting mesh can then be achieved using the Quadric Edge Collapse Decimation technique (Garland and Heckbert, 1997). Subsequently, the simplified and smoothed STL meshes were converted to CAD models using FreeCAD ⁵ . During this step, the common parts of the cell membrane between two intersecting cells are removed by the conversion algorithm, i.e., the cell cytoplasms are not separated by a membrane. The CAD geometries were, thereafter, utilised to generate quality tetrahedral volume meshes in Netgen.

Using 27 images from different cartilage zones, the classical segmentation approach together with the ellipsoidal fitting yielded an average cell volume of 2,157 ± 957 µm³. This cell volume is consistent with the reported value of 2,270 ± 590 µm³ by Lv et al. (2019), which we used as a benchmark. It indicates that the segmentation process is reliable. Furthermore, the obtained geometric model comprising the fitted ellipsoids resembled the original microscope image very well (see example in Figure 2).

However, some of the processed geometries could not be meshed. We identified small voids between intersecting or touching cells as the culprit. The error could be fixed by adding small spheres to cover the void. For that, an initial intersection check was performed on all objects and then an indexed list of the touching objects was created. The indices were mapped to a graph employing NetworkX (Hagberg et al., 2008) ⁷ . The nodes corresponding to intersecting objects were connected. In this way, cycles (loops) within the graph could be detected. Whenever a cycle in the graph is detected, the common volumes between the individual cells are computed and fused. In the case of a void, the fused common volumes do not form one single solid. The located void spaces were filled with a sphere.

The training of StarDist on our workstation without GPU support could not finish due to memory issues. Still, it succeeded on the HPC post-processing node after around five hours for 400 epochs when the loss value had converged. Subsequently, the trained deep learning model was validated with test images of varying quality. For that, we chose a good-quality image, two noisy images, and one highly noisy image. The predictions of the deep learning model, which was trained on our data, are shown in Figures 3, 4. The results suggest that the deep learning model performed well not only on high-quality images but also on noisy and highly noisy images. The precision, recall and F1 score with a recommended IoU threshold τ of 0.5 on the validation and test data, given in Supplementary Table S2, indicated that no over- or under-segmentation could be observed as those performance metrics are larger than 0.7. The evaluation metrics at different IoU threshold are provided in Supplementary Figure S4, S5. Moreover, the model achieved a test accuracy of approximately 80%, while the validation accuracy was around 60%. This result is in contrast to the general expectation that the performance metrics of the validation should be higher than those values from test data. We believe that our result can be explained by the data quality. In the cropped images, which we used for validation, several objects touched the image boundary. StarDist cannot detect objects in contact with the boundary of the images well, which could explain the relatively low validation accuracy. In practice, the deep learning model achieved a good prediction despite the relatively low validation accuracy.

To underline this point, we considered the total number of cells and the cell volume as additional model evaluation metrics. The predicted cell volumes did not deviate significantly from those obtained from the classical ellipsoid-fitting approach (see Table 1 for a comparison of the obtained volume ratios). Nonetheless, the number of cells counted using the deep learning method was 9.5% lower than in the classical approach. Moreover, a visual comparison of expected output versus deep learning prediction (refer to Supplementary Figure S6) also shows that most of the cells can be detected but some of them have slightly smaller volume than the cell in the ground truth. Using four test images, the classical approach detected 748 cells with an average volume of 1,714 ± 1,180 μm³, while the deep learning method identified 683 cells with an average volume of 1,601 ± 620 µm³. It is worth noting that the average cell volume is smaller than the reported value for all samples due to the predominant presence of cells from the top zone in the test image. Typically, chondrocytes in the top zone have the smallest volume compared to cells in other zones.

The ilastik model predicted 1,050 cells with an average volume of 1,487 ± 1,343 μm³, resulting in a 28% increase in the number of cells compared to the classical approach. As can be observed from Figure 5, utilising the ilastik model resulted in the detection of certain cell parts as separate objects, leading to an increased number of cells and smaller average cell volume. The StarDist prediction demonstrated better concordance with the classical approach, as evidenced by a smaller discrepancy in both the number of cells detected and their average volume. Moreover, due to the rugged shapes of the cells segmented using ilastik, it is necessary to apply post-processing techniques such as hole-filling, dilation, and erosion to obtain smooth and complete cell borders before converting them into STL geometry. Hence, StarDist is the method of choice in our specific case.

With CellPose3D, all four predicted images showed only a black background instead of cells (see Supplementary Figure S2). Due to the significant increase in training and prediction times compared to StarDist (approximately five times for training and eight times for prediction), we did not conduct further investigation into the algorithm failure in our case. It is worth noting that StarDist, unlike CellPose3D, does not require a data preparation step and uses the original format of the images as input. StarDist appears to be a more suitable algorithm for our specific application.

In order to better assess the performance of different segmentation approaches, we analyzed the distribution of cell volumes (Supplementary Figure S3). A few artifacts having a volume greater than 5,000 µm³ were identified as cells by the advanced classical approach. These artifacts shared the same pixel intensity as cells and could not be separated by means of the Distance Transform Watershed algorithms. Consequently, based on their fitted ellipsoid volume, these artifacts were excluded from the final geometry during the defining geometry input step. On the other hand, when utilising StarDist, all cell volumes fell within a reasonable range, specifically smaller than 4,000 µm³. It appears that StarDist can avoid wrong predictions of objects that deviate significantly in size from the training data. However, this characteristic may also be viewed as a disadvantage, because it necessitates ensuring that new images have a similar pixel size to the training images. Lastly, as previously mentioned, the usage of ilastik resulted in the detection of certain cell components as distinct objects. Hence, this led to an increased count of cells and a smaller average cell volume.

Regarding the model configuration of StarDist, we followed the advise by the developer. The training patch size was set to (16, 512, 512) to guarantee that it is smaller or equal to the size of the smallest annotated training image. Additionally, the grid parameter of (1,4,4) was used to specify the maximal size of well-segmented objects based on the neural network’s receptive field. The training batch size could not be set greater than one due to a lack of memory. A flexible parameter in the fitting procedure is the number of rays. We found that an adequately accurate reconstruction of the 3D cell in ground-truth images can be acquired with 96 rays. The evaluation metrics of the model did not improve when the number of rays was increased to 128. An even higher number of rays led to increased memory consumption that was beyond the capabilities of our systems. Moreover, we tuned the learning rate in conjunction with different model depths. The detailed results are provided in Supplementary Table S1. Based on our analysis, increasing the model depth to four resulted in memory issues. Additionally, when using a learning rate of 0.1 and 0.01, the model did not show any improvement as the loss fluctuated and remained relatively constant (refer to Supplementary Figure S7, S8). We observed that the model achieved the best performance with a learning rate of 0.0003 and a U-Net depth of two. As can be seen from Supplementary Table S1, reducing the learning rate to 0.0001 gained 1% in validation accuracy, but the training accuracy was reduced by 3%. Moreover, the model trained with a learning rate of 0.0003 achieved a lower training loss (0.5042) compared to the model trained with a learning rate of 0.0001 (0.5576) while the validation loss remained almost the same (0.6744 and 0.6713). Therefore, the learning rate of 0.0003 performed better in terms of training and is the value of choice. Our selection was further confirmed by the testing accuracy: 80% with the learning rate of 0.0003% and 77% with the learning rate of 0.0001.

By initialising the model with pre-trained weights from the StarDist model Rasse et al. (2020) along with the optimal hyperparameters mentioned earlier, all the accuracies of the models were reduced. Specifically, the training, validation, and test accuracies were approximately 2%, 3%, and 4% lower than the values when training the model from scratch. This decline in performance can be attributed to the significant influence of initial weights on the convergence of neural networks. Ill-suited weight initialisation might hinder the model’s ability to adapt and learn effectively from the new dataset (Glorot and Bengio, 2010; Sutskever et al., 2013).

We chose fitted ellipsoid images as the ground truth data because they were in good agreement with the manual segmentation (see Supplementary Figure S9) and StarDist is specialised for objects with blob-like characteristics. Alternatively, we used the images after applying the auto-threshold (see an example in Supplementary Figure S1) as ground truth. However, they produced poor predictions because the cells could have holes, rugged shapes or were incomplete close to the image boundary because several cells were not dyed completely. The training in StarDist involves fitting the ground-truth labels with star-convex polyhedra. Eventually, we concluded that the deep learning-based segmentation using StarDist yielded sufficiently good results to prepare geometries for the FEM simulations. Strikingly, the deep learning prediction using the trained model took less than three minutes on a standard workstation while the advanced classical technique, including the pre-processing and ellipsoidal fitting, took approximately 15 minutes to process one image. By utilising the HAUMEA HPC cluster’s post-processing node, the prediction time for the deep learning method was reduced to less than one minute, whereas the classical method required seven minutes to process one image.

The STL geometries generated from deep-learning-predicted images maintained the total volume and surface area with errors of less than 0.1% (Table 2). However, simplification is necessary to convert the generated STL meshes into CAD models. Table 1 provides the number of faces for the original STL geometry, as well as its simplified and smoothened version. The initial STL geometries were remeshed to generate an isotropic mesh. This process maintains the original volume and surface area. To ensure successful tetrahedral meshing, the obtained meshes were subjected to smoothening and subsequently simplified. Smoothing primarily contributed to a reduction in the surface area, while simplification was mainly responsible for a decrease in volume. As shown in Table 2, the final CAD models’ total cell area and volume decreased by approximately 15% and 5% with respect to those values from the deep learning predicted images. Smoothening the stepped structure of the STL geometry and the shape edges at the top and bottom of the cells due to the reconstruction of the marching cube algorithm from 3D stack images led to a large surface reduction. Nonetheless, the smoothening and simplification steps are inevitable to achieve FEM meshes.

We tested three different smoothing techniques in PyMeshLab, including Laplacian smoothing (surface preserving), Taubin smoothing, and two steps smoothing that included normal smoothing and vertex reposition. Among these techniques, the Taubin algorithm was the most suitable for preserving the cell volume. Crucially, it is the only technique that allows for successful FEM meshing. Moreover, two simplification methods available in PyMeshLab, edge collapse for marching cube meshes and quadric edge collapse decimation, were evaluated. The latter was found to be more effective in preserving the original surface area and volume of the STL mesh. The smoothened and simplified STL geometries using various smoothing techniques are illustrated in Figure 6.

The convergence of FEM meshes generated from deep-learning-based geometries was observed as the numerical simulation produced consistent results (with relative differences less than 0.05% as shown in Supplementary Figure S20) upon increasing the number of elements in the mesh through the use of smaller mesh sizes. The time required to convert the original STL generated from deep learning prediction to a CAD model is approximately five minutes while creating a CAD model from ellipsoid parameters derived from the classical segmentation approach takes less than one minute. Thus, using a deep learning approach on a standard workstation, creating a CAD model from the original image takes about eight minutes, whereas using the classical segmentation approach takes around 16 minutes.