Contributions of deep learning to automated numerical modelling of the interaction of electric fields and cartilage tissue based on 3D images

Electric fields find use in tissue engineering but also in sensor applications besides the broad classical application range. Accurate numerical models of electrical stimulation devices can pave the way for effective therapies in cartilage regeneration. To this end, the dielectric properties of the electrically stimulated tissue have to be known. However, knowledge of the dielectric properties is scarce. Electric field-based methods such as impedance spectroscopy enable determining the dielectric properties of tissue samples. To develop a detailed understanding of the interaction of the employed electric fields and the tissue, fine-grained numerical models based on tissue-specific 3D geometries are considered. A crucial ingredient in this approach is the automated generation of numerical models from biomedical images. In this work, we explore classical and artificial intelligence methods for volumetric image segmentation to generate model geometries. We find that deep learning, in particular the StarDist algorithm, permits fast and automatic model geometry and discretisation generation once a sufficient amount of training data is available. Our results suggest that already a small number of 3D images (23 images) is sufficient to achieve 80% accuracy on the test data. The proposed method enables the creation of high-quality meshes without the need for computer-aided design geometry post-processing. Particularly, the computational time for the geometrical model creation was reduced by half. Uncertainty quantification as well as a direct comparison between the deep learning and the classical approach reveal that the numerical results mainly depend on the cell volume. This result motivates further research into impedance sensors for tissue characterisation. The presented approach can significantly improve the accuracy and computational speed of image-based models of electrical stimulation for tissue engineering applications.

• A Gaussian filter was applied to a region of interest (ROI).
• A filter width was selected based on voxel sizes.
• A sphere diameter was specified based on the minimum expected object diameter in the dataset to perform background subtraction.
• If background subtraction is not selected, a absolute intensity threshold is applied.
• A cell volume was calculated after filling the holes and removing manually the outliers, such as small or large objects.

Trainable Weka Segmentation
Initially, we tried a trainable machine learning approach.We used the Trainable Weka Segmentation (TWS) method (Arganda-Carreras et al., 2017), which is well-established for microscopy pixel classification and available in the popular software ImageJ.First of all, because only one image can be loaded when training data, we chose ten slices from four images (two good-quality and two noisy images) to create one single image.This approach aimed to encompass a wide range of pixel intensities present in our dataset.A set of training pixels (STP) for each class (cell and background) were manually annotated on different slices by drawing traces using the Graphical User Interface (GUI).Although an acceptable classification with a small number of annotations can be obtained (Arganda-Carreras et al., 2017), we have marked 80 STPs for each class.In particular, the image has 80 STPs marked at different slices of 3D images using the cursor on the interactive GUI.In the setting panel, training features were kept as default in TWS, and image segmentation is obtained employing the Fast Random Forest algorithm (Breiman, 2001).The training was conducted on a workstation with Intel Core i7-10700 (2.90 GHz, 16 MB) and 64 GB of RAM.
The TWS training time is roughly 5 hours.The TWS classifier demonstrated excellent performance with high-quality images.However, for noisy images, the model produced relatively unsatisfactory results (as shown for example in Figure S1).Further, the TWS cannot separate touching objects or connect different parts belonging to one cell due to uneven distribution of dye, despite several attempts to manually mark the ground truth as precisely as possible.Consequently, the Distance Transform Watershed is required to apply once the image is fully segmented.Due to the satisfactory results achieved by the classical approach, we did not further fine-tune the TWS classifier.

DeepImageJ
An alternative method is to segment 3D fluorescent images using deep learning.Many software libraries and applications have surfaced in recent years.After an initial assessment, we focused on a selection of software solutions that in our opinion covered the field well enough.We assessed pre-trained deep learning We searched for 3D pre-trained models online2 .It turned out that most of the available models are either for electron microscopy (EM) or for detecting cell boundaries.We found only two models that were used for cell segmentation from confocal microscopy (Wolny et al., 2020;Bondarenko et al., 2022).The aforementioned deep learning models utilised a 3D U-Net (Ronneberger et al., 2015).Additionally, the three most popularly downloaded EM models were tested.However, neither of these models could be successfully applied to our data.We believe that the main reason for the failure is due to the differences between the training data and our data.The images used for pre-training were from different cell types in different tissues (such as cartilage, plants, and mouse embryos) and were acquired using different dyes.Thus, the training images may have had a different contrast, noise, and background intensity compared to our data.

Cloud-based Deep Learning Platforms
As the pre-trained models did not turn out to be capable of segmenting our data, we decided to train deep learning models on our specific dataset.
Biomedisa (Lösel et al., 2020) andZeroCostDL4Mic (von Chamier et al., 2021) were chosen because they are examples of cloud-based solutions and do not require own hardware.Biomedisa supports the training of CNN for the segmentation of 3D image data based on a 3D U-Net and is accessible through a web browser connected to a computing cluster.ZeroCostDL4Mic is a platform hosted on Google Colab and offers multiple segmentation algorithms, including also StarDist.A free version of Google Colab offers two Intel Xeon CPUs with one core and 12 GB RAM, and an NVIDIA Tesla K80 GPU with 12 GB VRAM and a runtime of up to 12 hours.A step-by-step instruction was outlined in the self-explanatory Jupyter Notebooks within the ZeroCostDL4Mic and in Biomedisa's website3 , eliminating the need for extensive computational expertise.
However, in both cases, we experienced out-of-memory errors indicating that the provided cloud hardware was not sufficient to train the provided algorithms with our data.

CellPose3D
Cellpose3D is an extension of the Cellpose algorithm (Stringer et al., 2020) designed for 2D application.It utilises the full 3D information considering the spatial context and volumetric characteristics of cells to enhance the segmentation accuracy on 3D image data.It also offers simplified training data preparation and instance reconstruction compared to the original version.We also investigated transfer learning on a pre-trained Cellpose3D model.The model was initially trained on 125 3D stacks of fluorescent confocal microscopy images depicting cell membranes in the meristem of the plant model organism A. thaliana.We fine-tuned the model with our dataset.The CellPose3D model underwent training and testing using identical training and testing data as the StarDist model.The data was prepared according to the instructions4 to meet the requirements of the CellPose3D model, which includes formatting and storing the data in an HDF5 file format.The default configuration settings were maintained throughout the training process, conducted on the post-processing node of the HAUMEA HPC cluster.The CellPose3D involves two primary steps in prediction.The network is first used to generate predictions for gradient fields and foreground segmentation based on the image data.Post-processing techniques are then employed to refine network predictions and reconstruct instance segmentation.The original code, training pipelines, and application pipelines are openly accessible at https://github.com/stegmaierj/Cellpose3D.Our modified scripts regarding to our specific application is given at dx.doi.org/10.5281/zenodo.7930409.
The CellPose3D models were trained for approximately 20 hours until the loss values reached convergence after 300 epochs.In all predicted images, no cells could be detected (see Figure S2).The failure of the predictions was due to the foreground-predicted images, as they identified almost the entire area as the foreground.Given the considerable increase in training and prediction times compared to StarDist as well as the additional data preparation step, we did not pursue further investigation into the algorithm failure.

ilastik
ilastik offers machine-learning-based bioimage analysis with a user-friendly interface.Its interactive workflow does not require a large 3D training dataset.Machine learning models for cell segmentation can be trained efficiently and intuitively in ilastik without extensive coding experience.The ilastik segmentation pipeline involved two main steps using its interactive GUI5 .First, we performed pixel classification by training a random forest classifier to distinguish between background and cytoplasm.We manually annotated voxels until additional training data no longer improved segmentation results on the training data.Training features were optimised using automatic feature suggestions provided by ilastik.Second, we identified individual cells by applying the hysteresis thresholding method to the results of the previous step.This step is an unsupervised technique that does not require any training data.Following the examination of various threshold value combinations, we determined that a threshold of 0.85 and 0.4 provided the most appropriate results for our specific case.The high threshold (core) and low threshold (final) are employed to differentiate and segment touching objects for detection purposes.The main paper includes a description of the results, providing a comparison with the StarDist algorithm.
Table S1.The effect of learning rate and the U-Net depth on the validation accuracy at the IoU threshold of 0.5.N.A indicates that a prediction process failed to due to an out of memory.

U-Net depth
Learning rate 0.001 0.0003 0.0001 0.00005  Using the advanced classical approach, a few artifacts with volumes exceeding 5000 µm 3 were misidentified as cells.These artifacts had similar pixel intensities to cells and couldn't be filtered based on volume alone.They were excluded from the final geometry using ellipsoid volume measurements.In contrast, StarDist produced reasonable cell volumes, all below 4000 µm 3 .This approach avoids predictions of objects that deviate significantly in size from the training data.ilastik detected certain cell parts as separate objects, leading to an increased cell count and smaller average cell volume.Table S3.Model parameters for UQ given as a function of the uniform distribution U .x, y, z: coordinates of the cell, V scale : scaling factor of the volume, θ: angle of the cell with respect to its local axis, dm: membrane thickness, σm: conductivity of the cell membrane, σ buf : conductivity of the cell culture medium, σcyt: conductivity of the cytoplasm, ε m r : relative permittivity of the cell membrane, ε buf r : relative permittivity of the cell culture medium, ε cyt r : relative permittivity of the cytoplasm

Case
Parameter Distribution Reasoning x / µm U(−10, 10) Assumptions for position error in x-axis y / µm U(−10, 10) Assumptions for position error in y-axis Case 1 z / µm U(−10, 10) Assumptions for position error in z-axis V scale U(0.8, 1.2) Scaling factor of the cell volume θ / °U(10, 90) Rotation angle to cover orientation error d m / nm U(3.5, 10.5) Assumptions from (Ermolina et al., 2000) σ m / µS m −1 U(0.08, 56) Assumptions from (Ermolina et al., 2000) σ buf / S m −1 U(0.95, 1.05) Assumptions from (Fear and Stuchly, 1998) Case 2 σ cyt / S m −1 U(0.033, 1.1) Assumptions from (Ermolina et al., 2000) ε m r U(1.4,16.8) Assumptions from (Ermolina et al., 2000) ε buf r U(60, 80) Assumptions from (Escobar et al., 2020) and (Zimmermann et al., 2022) Assumptions from (Ermolina et al., 2000) Figure S11.An illustration of geometries for UQ of the geometrical variation.One cell was always fixed at the top right corner of the box.A copy of the fixed cell was modified as shown in the three panels from left to right: the cell was placed in the centre of the square box; the cell was scaled to 80% of its volume, rotated by 45°and located at x = y = z = 10 µm; the cell volume was increased by 20%, rotated by 90°a nd located at x = y = z = −10 µm.The parameter combination in the above geometries was chosen for visualising examples of the geometrical models used in the UQ approach.The actual combination of cell location, rotation angle and volume scaling factor used in the UQ runs relies upon the algorithm and parameter distribution.

Figure S12
. The mean value and 90% prediction interval of the real part (left) and the imaginary part (right) of the impedance are demonstrated over a frequency range.The frequency-dependent first order Sobol indices are presented for parameters with a Sobol index of more than 0.05.The assumptions for the UQ analysis of the geometrical parameter, namely x, y, z position, volume and cell angle, are given in Table 1.Note that the other cellular dielectric properties were set at a constant.
)UHTXHQF\+] 5H=Ω 0HDQ SUHGLFWLRQLQWHUYDO )LUVWRUGHU6REROLQGLFHV σbuf εbuf Sobol indices are presented for parameters with a Sobol index of more than 0.05.The assumptions for the UQ analysis of the cellular dielectric parameter and the membrane thickness are given in Table 1.Note that the cell volume, angle and position were set at a constant.
. The mean value and 90% prediction interval of the conductivity (left) and the permittivity (right) are demonstrated over a frequency range.The frequency-dependent first order Sobol indices are presented for parameters with a Sobol index of more than 0.05.The assumptions for the UQ analysis of the cellular dielectric parameter and the membrane thickness are given in Table 1.Note that the cell volume, angle and position were set at a constant.

Figure S1 .
Figure S1.2D view of the middle slice from a noisy image.From left to right: the original image, the result of the TWS segmentation and the result of the classical segmentation after applying auto-threshold (the thresholding image), respectively.The TWS model performed poorly on noisy images (indicated by the blue ellipse) and was unable to separate touching objects or connect different parts of a cell (marked by the green ellipse), despite manual ground truth marking.As a result, the use of Distance Transform Watershed is necessary for post-segmentation processing.

Figure S2 .Figure S3 .
Figure S2.From top left to top right: the original image of A1S1, the final prediction from CellPose3D, and the foreground prediction from the neural network.From Bottom left to right: the neural network's predicted gradient maps in x,y and z-direction, respectively.The middle panel on the top row is the final predicted image reconstructed from the foreground prediction and the predicted gradient maps.No cells were detected in the final prediction.Its failure can be attributed to foreground identification, as a significant portion of the image area was erroneously classified as foreground (in white).

Figure S6 .
Figure S6.Superimposed positions of the expected output (the yellow ellipsoids) onto the deep learning created masks.The yellow ellipsoids highlight the regions that are expected to be segmented, while the deep learning created masks show the actual segmentation results.The majority of cells were detected, although some of their volumes were slightly smaller compared to the ground truth.

Figure S7 .
Figure S7.Training loss curve during the training and validation process with the learning rate of 0.1 and the U-Net depth of 2. The loss did not show any improvement as it increased after 50 epochs and remained relatively constant.

Figure S8 .
Figure S8.Training loss curve during the training and validation process with the learning rate of 0.01 and the U-Net depth of 2. The training loss is better than the case with the learning rate of 0.1, but still both training and validation loss did not exhibit any significantly improvement.

Figure S9 .
Figure S9.Overlay image of fitted ellipsoids and the original cells.The ellipsoids effectively represent cell structures.

Figure S10 .
Figure S10.An example of a slice in a training pair.Left: Input data as an original image.For visualization purposes, the Look-Up Tables (LUT) in ImageJ were modified to the "Fire" color scheme.Right: Ground true as a fitted ellipsoid image.

Figure S13 .
Figure S13.The mean value and 90% prediction interval of the real part (left) and the imaginary part (right) of the impedance are demonstrated over a frequency range.The frequency-dependent first order Sobol indices are presented for parameters with a Sobol index of more than 0.05.The assumptions for the UQ analysis of the cellular dielectric parameter and the membrane thickness are given in Table1.Note that the cell volume, angle and position were set at a constant.

Figure S15 .Figure S16 .
Figure S15.The mean value and 90% prediction interval of the real part (left) and the imaginary part (right) of the impedance are demonstrated over a frequency range.The frequency-dependent first order Sobol indices are presented for parameters with a Sobol index of more than 0.05.On the contrary to the results in Fig.S13, we neglected the uncertainty in the dielectric properties of the ECM and used fixed values.The assumptions for the UQ analysis of the cellular dielectric parameter and the membrane thickness are given in Table1.Note that the cell volume, angle and position were set at a constant.

Figure S17 .
Figure S17.The impedance plot of the A4S1 sample using the deep learning approach (dashed line and denoted by DL) and the classical approach (straight line).

Figure S18 .
Figure S18.The impedance plot of the A5S2 sample using the deep learning approach (dashed line and denoted by DL) and the classical approach (straight line).

Figure S19 .
Figure S19.The relative difference of the conductivity and permittivity of an individual sample using classical and deep learning segmentation in the meshing workflow.

Figure S20 .
Figure S20.The relative difference of the computed dielectric properties using moderate and fine mesh size.

Table S2 .
The precision, recall, F 1 score and accuracy with IoU threshold τ of 0.5 for validation and test dataset of the StarDist model.