Following the method previously described, a transgenic Hydra vulgaris (strain AEP) line was established to express GCaMP7s in the cytosol and tdTomato in the nucleus of cells derived from the interstitial stem cell lineage. The plasmid (Addgene catalog no. 102558) was modified by replacing the actin promoter with the EF1-alpha promoter and inserting DsRed downstream of the GCaMP6s sequence. To ensure the nuclear localization of DsRed, a P2A self-cleaving peptide sequence containing a nuclear localizing signal (cccaagaagaagaggaaggtg) was inserted between GCaMP6s and DsRed. The nuclear localization of DsRed was confirmed by electroplating the plasmid into Hydra. To enhance fluorescence intensity, GCaMP6s was replaced with jGCaMP7s and DsRed was replaced with tdTomato. Finally, for microinjections, the EF1-alpha promoter was replaced with the Actin promoter. All fluorescent reporter gene sequences were codon optimized specifically for Hydra. Standard embryo microinjection of this plasmid was performed and transgenic hatchlings expressing cytoplasmic GCaMP7s and nuclear tdTomato in the interstitial cell lineage were isolated. Transgenic animals were bred until all neurons were expressing both cytoplasmic GCaMP7s and nuclear tdTomato. Transgenic Hydra were cultured using standard methods in Hydra medium (1 mM calcium chloride dehydrate, 0.33 mM magnesium sulfate anhydrous, 0.5 mM sodium bicarbonate, 0.03 mM potassium chloride) at 18 °C on a 12 h light/12 h dark cycle. They were fed freshly hatched Artemia nauplii twice per week.

 

 

Dual-labeled transgenic Hydra were prepared for imaging as described. Imaging was performed using a custom dual-channel spinning disc confocal microscope (Solamere Yokogawa CSU-X1) with a sCMOS camera for each channel (Teledyne-Photometrics Prime-BSI). Samples were simultaneously illuminated with both 488 nm and 561 nm lasers (Coherent OBIS). Emission light was split with a dichroic mirror sending green light to one camera and red light to the other. Cameras were aligned with pinholes and images were registered during processing. Images were captured with a frame rate of 10 frames per second using a 6X objective (Navitar HR Plan Apo 6X/0.3) and Micro-Manager software.

 

 

Immunostaining was performed using the previously established protocol optimized for double-labeling experiments in Hydra with the anti-Hydra cadherin antibody, using the following antibodies: Primary antibodies: Hydra cadherin antibody (rabbit) 1:1000 (Thomas Holstein, Heidelberg), anti-tdTomato (goat) 1:200 (Origene, AB8181-200); Secondary antibodies: Alexa 488 donkey anti-rabbit 1:1000 (ThermoFisher A-21206), Rhodamine Red-X donkey anti-goat 1:250 (Jackson 705–295-147). Fixed and stained animals were imaged on the same spinning disc confocal setup described above. Image stacks were taken with a 40X objective (Olympus LUMPlanFl/IR 40X/0.8 W) every 0.5 μm and processed with ImageJ. Images shown are maximum intensity projections of short stacks (5–20 μm).

 

 

 

To learn a segmentation model and validate our experiments, we sampled 20 tdTomato images of Hydra on various positions of the animal (e.g., contracted, elongated). We manually annotated these images using the ImageJ draw tool resulting in 12,705 segmented neurons (~ 600/image).

 

To measure the segmentation performance, 5 random images were kept as a testing dataset. To account for performance variability, we repeated this process 5 times with different seeds, resulting in different training and testing images. We measured the performance with each seed and report the results as (mean ± std).

 

We used standard instance segmentation metrics (precision, recall, f1-score). But we did not base these metrics on IOU association, where a predicted instance is considered to be a true positive if it overlaps with a ground truth on a sufficiently large area. We rather based these metrics on distance association, where a predicted instance (spot) is considered valid if its mass-center is sufficiently close to the mass-center of a ground truth instance. This has two benefits: first, it is more meaningful for tracking as we only use the position of each detected instance. Second, it allows fair comparison with detection methods that do not try to produce a faithful segmentation or bounding box of the instances like the wavelet thresholding method.

 

The remaining 15 images were used to train and tune all the parameters of the segmentation model. We also measured the impact of the training data quantity by training the model with fewer images, from 1 to 15.

 

We compared two detection methods: an analytical method based on wavelet decomposition and thresholding, and a deep learning approach, StarDist.

 

 

The method has only two hyperparameters: the scale of the spots and the noise threshold. Both can be easily tuned by hand on a single image. Nonetheless, to be fully autonomous, the system performs a grid search on the training images and chooses the best performing ones. No further training is required. Moreover, to reduce the number of false positives, the neurons with less than 5 pixels area were filtered out.

 

 

The deep neural network must be trained and the hyperparameters tuned. We therefore split the training images into training and validation as most standard deep learning approaches do. 20% of the training images were used as validation images. With less than 5 training images, one image was kept for validation. Therefore, the procedure that we show here requires at least 2 training images, one to train the weights of the network, the other to validate the hyper parameters.

 

The official implementation of StarDist (https://github.com/stardist/stardist) was used to train, validate, and evaluate the performances with our dataset. For tracking, we provide a wrapper to perform the detection process using a trained StarDist model (https://partage.imt.fr/index.php/s/npwHJHZebxqGMPi).

 

 

The eMHT algorithm is implemented on the open-source imaging platform Icy (plugin Spot Tracking), which is coded in Java language. Instead of translating all the code into Python, we decided to directly call the Spot Tracking plugin in Icy in a headless mode from the TraSE-IN platform. To ensure compatibility with Icy, we implemented different functions for inputs (detections) and outputs (tracks) wrapping.

 

We believe that calling implemented tracking solutions in major imaging platforms such as Icy and ImageJ is an efficient and sustainable solution because the main imaging platforms are open-source, the tracking plugins (e.g., Spot Tracking in Icy, Trackmate in ImageJ) are regularly updated by developers, wrapping input and output variables in TraSE-IN is much less tedious and prone to implementation errors than translating all the Java code of tracking methods in Python, and TraSE-IN will easily integrate new tracking solutions developed on other platforms, if needed.

 

 

To stitch tracklets obtained with eMHT tracking, we implemented a cost-based algorithm. Briefly, a cost between all the tracklets is computed, and stitched tracklets correspond to the minimal global cost of all tracklets’ linking. Computing the minimum global cost corresponds to a linear assignment problem that we solved using the Jonker-Volgenant algorithm. The cost between tracklets was computed following the EMC algorithm. Input parameters of the stitching function (byotrack.implementation.refiner.stitching.emc2) are the smoothness α of the Thin Plate Spline interpolation and the non-linking cost η (paid for unlinked tracklets). In our tracking scenario (Fig. 4C), we set α = 10 to provide enough regularization to be robust to eMHT linking errors and set non-linking cost η = 5 pixels.

 

 

To solve the misalignment between nuclei and calcium signal, we implemented a sub-ROI tracking method. For each nucleus, we extracted a 25 × 25 pixels ROI centered on each tracked nucleus. Then, to find neuronal candidates in the GCaMP7s channel, we iteratively computed η local maxima within each ROI (η = 5 typically) and identified the most probable calcium signal. For this, we first applied a gaussian smoothing (σ = 1 pixel) and first identified the global maximum in the ROI. We then deleted the pixels in the neighborhoods of this maximum and reiterated the process η times. To bias extracted maxima towards the nucleus position (and avoid selecting an outlier GCaMP7s signal) we weighted the ROI intensities with Gaussian prior centered on the nucleus position (σ = 5 pixels).

 

In the first frame, the selected position of the calcium spot is the global maximum of the ROI (after having applied a Gaussian prior to bias intensities towards the nucleus position). For the following frames, we selected the closest maximum to the previous calcium position and updated the calcium position with an exponential moving average to attenuate the potential consequence of a wrong association. Finally, we extracted the calcium intensity for each neuron as the mean intensity within the 5 pixels radius circle centered on the tracked calcium position.

 

 

 

We compared extracted calcium signals with a normal Gaussian distribution using the scipy.stats.normaltest method. Any signal whose normality can be rejected with p-value lower than a user-defined threshold was kept. The other signals that might correspond to non-neuronal cells are displayed in an interactive feedback step to allow users to decide which cells they want to keep or discard from the analysis. This feedback step allows the user to visualize each cell corresponding to a Gaussian fluorescent signal. Then, by observing the cell’s location, size, shape, and activity over the timelapse movie, users can decide if they would like to discard the cell or keep it for further downstream analysis.

 

 

We used the tdTomato intensity (extracted as the mean intensity in a 5 pixels radius circle centered on the nucleus position in the red channel) as a control signal. ICA between the control signal and the raw calcium signal was used to remove motion artefacts. ICA relies on scikit-learn FastICA method (we used a tolerance of 0.5). As the ICA algorithm is stochastic, we ran it several times and selected the best outcome as the one with the minimum correlation with the control signal.

 

 

To remove the remaining baseline in the calcium signal, we applied frequency filtering using a 5th order Butterworth filter and a critical period of 100 frames. In addition, we smoothed the signal using a rolling average of size 5. Both the average size and critical period can be adapted with user feedback.

 

 

We used the constrained FOOPSI implementation in the CaImAn Python library (OASIS method of order 2). In Hydra, calcium spikes correspond to single action potentials. Therefore, we added a clustering step to aggregate the action potentials (electrical spikes) predicted with the previous deconvolution FOOPSI algorithm (less than 5 frames away). The clustering consists of first temporal Gaussian blurring (σ=5 frames) of the estimated electrical spikes, followed by local maxima extraction. Finally, we filtered out electrical spikes with low probability/intensity using a user-defined threshold to reduce the number of erroneous predictions. Since these computed probabilities/intensities vary wildly between different signals, using a fixed threshold would not work across different signals. To rectify this, an adaptive threshold defined as the maximum probability minus 2 standard deviations was used. This method provided robust results across all cells, allowing spike extraction to be applied simultaneously across an entire dataset with minimal user input.