2.1.

Clinical Data and Equipment

This pilot clinical study was conducted between September 2020 and August 2021, with 18 patients with histological diagnosis of PCa, enrolled on a prospective clinical trial approved by the Research Ethics Board (NCT03378856).

Planning mpMRI (3D T2-weighted FSE, b2000 DWI, $+/-\mathrm{DCE}$) were obtained on a 1.5T Siemens Aera Magnetom (Siemens Healthineers, Erlangen, Germany) using surface coils. Voxels on acquired T2 images were $1\times 1\times 1{\mathrm{mm}}^3$, b2000 images (diffusion-weighted images $b$-value of $2000{\mathrm{s}/\mathrm{mm}}^2$) consisted of $2.6\times 2.6\times 5{\mathrm{mm}}^3$ voxels, and voxels on calculated ADC maps were $1.8\times 1.8\times 4{\mathrm{mm}}^3$. Positron emission tomography/computed tomography (targeting prostate-specific membrane antigen: PSMA-PET/CT, 18F-DFCPyL)²² images were also acquired in a subset. For intraoperative imaging, a bk3000 ultrasound system was used with a BK endocavity biplane transducer (BK Ultrasound, Herlev, Denmark). Brachytherapy procedures were assisted by a prototype interventional system (Invivo/UroNav, Philips Disease Management Solutions, Gainesville) in the early phase of clinical deployment.

We used a custom system consisting of navigation and optical components (Fig. 1). The optical component contains a dual source (671 and 785 nm, Semrock, New York), a spectrometer (EmVision LLC, Fl), and the custom EM tracked optical probe (EmVision LLC, Fl), designed to perform in-situ, minimally invasive characterization.⁴⁵ The use of this subsystem, controlled by customized MATLAB R2017b (Mathworks, Massachusetts) software, allows one to stimulate the tissue and detect energetical shifts due to inelastic scattering from light–tissue interactions, which is correlated to molecular vibration modes.¹⁹

Fig. 1

Clinical setup: (a) ultrasound system; (b) near-infrared laser and spectrometer; (c) EM tracking system; (d) EM field generator; (e) 3D Slicer navigation system; (f) closeup of the RS probe fiber bundle next to cannula; (g) Raman excitation fiber; (h) Raman detection fibers; (i) EM sensor; and (k) schematic illustration of optical setup.^20,45

The navigation component is connected to the ultrasound system and uses MRI-TRUS fusion to project structures segmented on MRI over the real-time TRUS images. It is also connected to the EM tracking platform (Aurora NDI, Waterloo, Canada) consisting of a control unit, an EM field generator placed over the patient’s pelvic region, and three six-degrees-of-freedom EM sensors: the first one placed on the template as a fix reference, the second fixed to the TRUS-probe holder to track the field of view of the current 2D image, and the third one, small enough to be integrated to the custom optical probe lumen, to reconstruct the probe while navigating to pre-identified targets. A 3D Slicer module was created for visualization and control of the subsystems.⁴⁶ More details on the system could be found in previous works.^20,47

2.2.

Workflow and Data Acquisition

HDR brachytherapy procedures include the following steps: (i) importing the preoperative images and predetermined contours into the interventional system, (ii) performing a 3D reconstruction from 2D TRUS image and segmenting the prostate, (iii) registering images and propagating the MRI contours to TRUS, (iv) inserting and reconstructing the catheters, and calculating the dose plan, and finally, (v) delivering the dose. All these steps were performed with the patient under general anesthesia and in lithotomy position. Details of the intervention (brachytherapy workflow), image registration, and EM tracking can be found in Ref. 48.

The optical acquisitions were performed during HDR brachytherapy procedures after the elastic image registration (defined below); once the optical measurements were completed, the brachytherapy procedure resumed. The general diagram of the connections and the flow diagram can be seen in Fig. 2.

Fig. 2

Flowchart and connections diagram of the systems involved in the intervention. The larger boxes divide the pre-, intra-, and post-operative phases of the procedure. The left column in the intraoperative section contains the equipment and activities inherent to brachytherapy. The right column contains the steps and systems added to the interventions to carry out the present study.

Gross tumor volumes (GTV) were segmented on T2 MRI, based on mpMRI and PSMA-PET images; prostate and urethra were segmented in preoperative images (using Eclipse, Varian Medical Systems, Palo Alto) and intraoperative images (using the interventional system). Additionally, T2 images and contours were elastically registered to the intraoperative TRUS, using a surface-based algorithm integrated into the prototype interventional system; for this process, the prostate contours (from MRI and TRUS) were automatically centered, then manually aligned following the urethra angle, and then automatically deformed looking for an optimal surface correspondence. TRUS images and propagated contours were exported to our custom navigation system, enabling visualization of deformed structures (e.g., GTV) projected over real-time TRUS.

Targets (e.g., GTV center of mass, healthy tissue far from the GTV) were pre-identified based on the preoperative MRI, and the optical probe was navigated to them supported by a coaxial needle (cannula). Once at the site of interest, the dual optical source stimulated the tissue, and the spectrometer captured the response signal (from 50 to 100 RS spectra per site). The coordinates of the inspected sites in the TRUS reference system were recorded, and a confirmation biopsy was taken at the same location.

Biopsy cores were fixed and processed according to standard histopathologic procedures observed by an expert to identify patterns on the stained sample slides.⁴⁹ According to the predominant and secondary patterns, a report was generated presenting the Gleason score (GS), the grade group according to the International Society of Urological Pathology (ISUP GG), and the high-grade tumor percentage (HG) for each biopsy core.⁵⁰ For instance, a slide presenting tumor tissue, consisting of 80% pattern 4 and 20% pattern 3, will be reported as GS: $4+3=7$, ISUP GG: $3/5$, and HG:80%.⁵¹ Due to the impact over treatment planning, classification algorithms for PCa are usually trained for detecting tumors with $\mathrm{ISUP}\mathrm{GG}>1$, so initially, we set our “$\mathrm{ISUP}\mathrm{GG}>1$” prediction task by labeling the observations with $\mathrm{ISUP}\mathrm{GG}\le 1$ (including $\mathrm{ISUP}\mathrm{GG}=1$ and benign tissue) as false and sites with $\mathrm{ISUP}\mathrm{GG}>1$ as true.^1,41,52

The RS signal processing consisted of averaging all the spectra from one site, removing autofluorescence and cosmic rays, standard normal variate normalization, and finally assigning each pixel of the spectrometer to a Raman shift (Fig. 3).⁴⁵ We obtained a single spectrum per site at the end of this process. This was applied for the fingerprint (FP) and the high wavenumber (HW) region of the RS, according to the source used to stimulate the tissue. For this study, every single Raman shift on the spectrum was a feature for the next steps.⁴⁷

Fig. 3

Sample signal processing. (a) About 50 raw fingerprint acquisitions from one inspected site and (b) the resulting processed Raman spectrum.

To extract MRI-based features (Rad), the PyRadiomics platform was used,⁴² which offers up to 120 features from different categories; based on literature, we selected 8 first order and 8 GLCM features (Table 1) to be calculated on T2, ADC, and b2000 mpMRI.^8,40,41,52 Using the transformation matrix employed to register MRI and TRUS, we applied the inverse process to find the corresponding coordinates on the mpMRI for each inspected site. We calculated the radiomics features in a 5-mm radius spherical volume at each position.

Table 1

Radiomics features extracted from the three mpMRI sequences, and their identifiers.

2.3.

Classification Model Training

Given the low number of patients, which is inherent to a pilot clinical study, we decided to limit the complexity of the classification model as a strategy to avoid overfitting issues. Thus, we trained SVM models for binary classification, with a linear kernel and a cost matrix that doubles the penalty to false negatives ($C=[0,1;2,0]$), using MATLAB R2017b. For training and validation, we followed a leave-one-patient-out cross-validation (LOPOCV) scheme (i.e., models were trained with data from all patients but one, which was used for validation).

In a high-dimensional dataset, especially with a limited number of observations and patients, feature selection is crucial. To perform this step, we set a maximum number of features to be selected ($\max \_\mathrm{\text{nf}}$) and applied a three-step selection approach:

This three-step process was applied after dividing training/validation sets, i.e., the observations of the patient left out were not considered during the feature selection. As a result, the number of features selected for each fold may vary slightly, so the value reported in the results is the average of all folds.

We describe below three sets of experiments focused on the classification potential of collected features (RS and Rad) to train the predictive model.

3.1.

Clinical Data

In total, the dataset consisted of 47 inspected sites, with the corresponding histopathological report for ground truth, and 4650 features (Table 2). Sample images of two histopathological slides are presented in Fig. 4.

Table 2

Clinical data.

Fig. 4

Sample image registration and co-location results. (a) GTV originally segmented on preoperative mpMRI is (b) projected over the real-time TRUS using a surface-based elastic registration algorithm; the inverse process allows identifying the MRI coordinates corresponding to the inspected sites. (c) A 3D model of the plan can be rendered. (d) Pathological image samples of benign prostatic parenchyma and (e) acinar adenocarcinoma GS: $3+4=7$, $\mathrm{ISUP}\mathrm{GG}=2$, and HG: 40% to 50%.

Elastic registration allowed the projection of the deformed GTV contours on real-time TRUS, which was used to set up to five sites of interest and for navigation support during the intervention; the inverse process allowed the co-location of the inspected sites on the mpMRI, post-intervention [Figs. 4(a)–4(b)]. According to the voxel size of mpMRI and the volume defined for radiomics extraction (blue sphere on Fig. 4), ranges of 428–482, 34–38, and 24–32 voxels were used from T2, ADC, and b2000, respectively.

3.2.

Feature Combination Experiment

Following a dichotomization of the collected data during the brachytherapy procedures ($\mathrm{ISUP}\mathrm{GG}>1$ criteria), these experiments consisted of $21\mathrm{ISUP}\mathrm{GG}>1$ and $26\mathrm{ISUP}\mathrm{GG}\le 1$ samples. The AUC results for the different combinations of features are shown in Table 3, and the ROC curves for the groups of features involved in the model with the best performance are presented in Fig. 5.

Table 3

Classification performance of SVM models (LOPOCV) for discriminating ISUP GG>1 and ISUP GG≤1.

Fig. 5

ROC curve and optimal point for discriminating $\mathrm{ISUP}\mathrm{GG}>1$ and $\mathrm{ISUP}\mathrm{GG}\le 1$, for models including feature selection ($\max \_\mathrm{\text{nf}}=10$). AUC: Area under the curve.

The $\mathrm{FP}+\mathrm{Rad}$ model showed the best overall performance, based on AUC, followed closely by three other combinations: $\mathrm{FP}+\mathrm{HW}$, rad, and $\mathrm{FP}+\mathrm{HW}+\mathrm{rad}$. From these four models, $\mathrm{FP}+\mathrm{Rad}$ and $\mathrm{FP}+\mathrm{HW}+\mathrm{Rad}$ yielded a low mean support vectors ratio (number of support vectors needed to train the model, divided by the total number of observations) compared to the other two models.

Prediction accuracy, sensitivity, and specificity were calculated at the optimal point on the ROC curve, which is the closest point to the (0,1) corner of the plot. These metrics also showed the superior performance of $\mathrm{FP}+\mathrm{Rad}$ (Table 3). Regarding the prediction accuracy, $\mathrm{FP}+\mathrm{Rad}$ outperformed all other models but one for more than 9%; $\mathrm{FP}+\mathrm{HW}+\mathrm{Rad}$, which also combined optical and image-based features, was the only model close to this accuracy (only 2% below) but also 10% below in sensitivity.

3.3.

Number of Features Experiment

Using the combination of features with the best overall performance ($\mathrm{FP}+\mathrm{Rad}$) identified in the previous experiment (Sec. 3.2), we trained models varying $\max \_\mathrm{\text{nf}}$ for the feature selection process. As explained in Sec. 2.3, the selected features for each fold (LOPOCV) could vary slightly, up to the set limit ($\max \_\mathrm{\text{nf}}$), so the actual number of features reported is the mean number of features selected for all the folds. According to AUC results and metrics calculated at the optimal point (Fig. 6), there is an optimal region between 6.1 and 8.3 used features, resulting when setting the $\max \_\mathrm{\text{nf}}$ at 7 and 10, respectively. Limiting the number of features used to train the classification model aimed to reduce overfitting given the small datasets and improve generalization capabilities; but the remarkable drop in specificity observed in the graph for more than nine features emphasizes the importance of the feature selection process.

Fig. 6

Classification performance of the $\mathrm{FP}+\mathrm{Rad}$ model in function of the number of features selected.

The features selected on each fold for the model that used on average 8.3 features were not always the same, so Figs. 7(b)–7(d) shows how many times (# of folds) each feature was selected. FP Raman spectra were averaged by class for visualization [Fig. 7(a)], and from those 1801 features (Raman shifts), five features were selected more than 10 times (i.e., during the feature selection applied over each fold): 994, 1007, 1334, 1766, and $1772{\mathrm{cm}}^{-1}$. Similarly, radiomics features were averaged by class [Fig. 7(c)], and two features were selected more than 10 times: r17 (Energy from ADC) and r46 (DifferenceEntropy from B2000), both selected on all the folds (18 times). This total of seven most frequently selected features could make a set of features with great potential for further studies.

Fig. 7

Average and standard deviation (a) of processed FP Raman spectra and (b) histograms of selected features for RS; average and standard deviation (c) of mpMRI radiomics and (b) histograms of selected radiomics features. The radiomics features corresponding to each identifier are given in Table 1.

It can be observed, as expected, that the valleys in the Raman spectrum were not selected, except for some features just at the base of a significant peak, which contained information on the width of the peak (e.g., 1108 or $1575{\mathrm{cm}}^{-1}$). As for the radiomics features, it can be seen that, apart from the two most frequently selected features, only a few were selected at least once (r18, r20, r39, and r44); according to the selection method, none of the T2 features, GLCM from ADC, and first order from b2000 contributed to the classification task.

3.4.

Prediction Task Experiment

We can observe that the different criteria to set the classes for the prediction tasks resulted in different numbers of labels and observations (Table 4). The three observations with $\mathrm{ISUP}\mathrm{GG}=1$ (Table 2) were $\mathrm{HG}<20\%$, so the last set of labels (high grade) did not include these samples.

Table 4

Classification performance using the most frequently selected features, based on different prediction tasks.

The models trained for each prediction task were compared in terms of general performance (ROC) and the performance at the optima point, and the results are shown in Fig. 8 and Table 4.

Fig. 8

ROC curve and optimal point for models trained with the seven most frequently selected features, for the three different prediction tasks ($\mathrm{ISUP}\mathrm{GG}>1$, $\mathrm{ISUP}\mathrm{GG}\ge 1$, and high grade).