The study was a physician-led safety trial with a double-blind design conducted at the International Neurological Restoration Center (CIREN in Spanish), Cuba, collaborating with the Center for Molecular Immunology (CIM in Spanish). The study was registered at clinicaltrials.gov/ct2/show/NCT04110678?term=neuroEPO&cond=Parkinson&cntry=CU& draw=2&rank=1 with the number NCT04110678.

Participants were recruited in the outpatient Movement Disorders Clinic at CIREN from February 2015 to July 2016. Twenty-six patients with a clinical diagnosis of PD were selected from 46. The inclusion criteria for the selected sample met the criteria of the UK Brain Bank, severity stages of Hoehn and Yahr (1967) I–II, age between 40 and 70, good response to dopaminergic medication: 30% of change of motor symptoms as rated on the Movement Disorder Society-Unified Parkinson Disease Rating Scale (Goetz et al., 2007), good general health, without depression or cognitive impairment. The exclusion criteria were stringent and eliminated patients with cognitive deterioration, psychiatric disorders, pregnant women, hypertension or any non-compensated diseases, sepsis, or treatment with other drugs. The detailed list of inclusion criteria can be found in Garcia Llano et al. (2021). The severity of the PD was evaluated by the Hoehn and Yahr scale.

The 26 patients were randomly divided into two groups: One group consisted of 15 patients who received a 1 mL dose of intranasal Neuroepo per week for 5 weeks. Another group included 11 patients who took a placebo.

The evaluation of the motor function was performed by a neurologist of the Movement Disorders Clinic at CIREN, using the Movement Disorder Society motor scale (Unified Parkinson’s Disease Rating Scale section motor; UPDRS III).

We applied the levodopa equivalent dose (LED) method to standardize the intake of pharmaceutical agents by the patients. The LED was calculated using the conversion proposed by Schade et al. (2020). See Table 1 for the demographic and clinical description of this sample.

The EEG, neuropsychological, and motor evaluations were performed at two points: at the baseline and 6 months later. The patients were studied in the “ON” state for all the assessments planned during the early mornings. The patients were advised to take breakfast and their regular medication before attending the evaluations.

Other assessments included clinical blood tests and blood pressure to follow the health status of the patients (data not shown).

The project had ethics approval from the institutional review board at CIREN. Written informed consent was obtained from all patients and caregivers following the ethics standards at CIREN and CIM.

Cognitive performance was assessed with a comprehensive neuropsychological battery for global cognitive screening, comprising 32 variables from 9 tests. The tests included the Mini-Mental State Examination (Folstein et al., 1983) and the Dementia Rating Scale (Mattis, 1988); the Rey Auditory Verbal Learning test (Bean, 2011), the subtest letter-number sequencing of the Working Memory Index of WAIS III (Wechsler, 1997), the Rey Complex Figure, copy and delayed recall (Rey, 1941), Delis-Kaplan verbal fluency (Delis et al., 2001), Trail-Making (Reynolds, 2002) the Stroop color-word Interference test (Rivera et al., 2015), and the Frontal Assessment Battery (Dubois et al., 2000).

The initial 32 variables were screened to leave only those with a significant difference in t-tests between the two measurement points for the Neuroepo and control groups. FDR set the threshold with both p and q levels set to 0.05. Consequently, the final six neuropsychological variables included are described in Table 2.

The mood of the patients was assessed by the Hospital Anxiety and Depression Scale (HADS) validated for the Spanish population (Herrero et al., 2003) only at the baseline to decide the final inclusion of the patients.

Eyes closed resting monopolar EEG was recorded from Fp1, Fp2, F3, F4, C3, C4, P3, P4, O1, O2, F7, F8, T3, T4, T5, T6, Fz, Cz, and Pz using the International 10/20 system, all referenced to the linked earlobes [A1 (left)-A2 (right)]. The EEG was recorded using a Cuban neurometric system Medicid IV system and the Track Walker TM v5.0 software¹ at the CIREN Neurophysiology Lab. For all electrodes, the impedance did not exceed 5 kΩ. EEG signals were amplified 10,600-fold, bandpass filtered from 0.5 to 30 Hz, and sampled by a 12-bit analog-to-digital converter at 200 Hz.

The EEG was recorded at two assessment points: at baseline before the intervention and 6 months after administering the drug or placebo. We selected 2.56-s artifact-free segments selected by an expert (TVA) for the analysis using visual inspection. For the removal of other artifacts due to interfering physical or physiological activities, we employed different semi-automatic plugins based on EEGLAB (Delorme and Makeig, 2004) under the Matlab platform. Because of artifacts and poor quality EEG recordings, one participant from the placebo group was eliminated from the final sample. Thus, the final sample for the analysis was Neuroepo n = 15, placebo = 10, total n = 25, reduced from the original n = 26.

The topographic (qEEG) and tomographic (qEEGt) quantitative EEG methods have been explained elsewhere (Taboada-Crispi et al., 2018; Bringas Vega et al., 2019; Bosch-Bayard et al., 2020) and are now implemented as the qEEGt plugin at the CBRAIN (Bosch-Bayard et al., 2020) platform, available from the CBRAIN portal: https://portal.cbrain.mcgill.ca.

For the Quantitative EEG analysis (qEEG), the EEG signals in time were re-reference to the Average Reference montage and corrected by the Global Scale Factor (Hernandez et al., 1994). This factor accounts for variability in the EEG due to factors unrelated to the neurophysiology, such as the skull thickness, hair volume, impedance, recording conditions, and others, which may introduce a baseline bias affecting the spectra from different participants. The EEG was then transformed to the frequency domain using the Fast Fourier Transform (FFT). An average of 24 epochs of EEG artifact-free and quasi-stationary signals, manually selected by the expert neurophysiologists, were used to calculate the Power Spectral Density (PSD), both at the scalp and well at the sources. Each epoch comprised 512-time samples (i.e., 2.56 s). The transformation to the frequency domain was carried out using Bartlett’s method (Møller, 1986) by averaging the cross-periodograms. This procedure yielded 49 cross-spectral matrices for the 19 electrodes in a frequency range from 0.78 to 19.14 Hz, with frequency bins every 0.39 Hz.

For the calculation of the PSD at the EEG sources, the VARETA (Variable Resolution Electromagnetic Tomography) methodology, as implemented in the Neuronic software (Bosch-Bayard et al., 2001) and CBRAIN plugin (Bosch-Bayard et al., 2020) was used. VARETA is a discrete spline EEG inverse solution based on a forward model that incorporates anatomical constraints using the template of ICBM (probabilistic brain atlas) created by the Montreal Neurological Institute (MNI) (Evans et al., 1993). Source spectra were then Log-transformed to achieve an approximately Gaussian distribution.

To recap, the data analyzed comprised these measures:

To deal with the high dimensionality of data and small sample size, we employed probabilistic Canonical Correlation Analysis (pCCA) to obtain cognitive and qEEG cognitive latent variables reflecting common sources of variation. We then implemented a simple causal (counterfactual) mediation model with repeated measures to investigate the dose-effect on qEEG and cognitive variables using these latent variables.

Whitening, or sphering, is a preprocessing step for data analysis and machine learning models. A set of random variables with a known covariance matrix is transformed into a latent space with identity covariance (Zhang, 1998). Whitening is used to achieve unique features in data and remove redundant information. It also fulfills the purpose of data compression (Kessy et al., 2018). Whitening can be formally expressed as:

Z is the observed variable, Z is the latent whitened variable, and Wz is the whitening matrix or unmixing matrix. A recent study implemented pCCA as a whitening transform to integrate high-dimensional gene expression and methylation data for lung carcinoma (Jendoubi and Strimmer, 2019). The purpose of using CCA for constructing Wz is twofold. It simultaneously whitens two variables while constraining the cross-correlation matrix to achieve highly correlated unique features between two measured variables. Jendoubi and Strimmer (2019) also demonstrated that their implementation of CCA can identify both positive and negative associations (unlike classic CCA, which allows only positive correlations).

We have implemented data whitening with pCCA to integrate qEEG and Cognitive variables using the whitening package in R (Strimmer et al., 2021). The first factor from whitened qEEG and Cognition (qEEG,Cog) was used for the mediation analysis.

To investigate if the effect of Neuroepo dosage on Cognition is mediated by qEEG, we implemented the mediation model shown in Figure 1.

In this model Dosage of Neuroepo is the treatment variable. The latent variables qEEG (qEEG) and Cognition (Cog) are the mediator and outcome, respectively. The significance of the paths is analyzed using the following general repeated-measures, mixed-effect models:

Here Y is the outcome variable, X and S are the fixed and random effects covariates. β and θ are estimated parameters for fixed and random effects. This general model was instantiated in two specific linear mixed effect models:

These models were fitted using the lme4 package (Bates et al., 2015) from R. We also performed a power analysis using the “simr” package in R, which runs the simulation-based power analysis for linear mixed effect models (Green and MacLeod, 2016). Power analysis for both models was carried out with 100 simulations. A statistical power value threshold of 80% or more was chosen to identify adequate statistical power.

Furthermore, we combined the estimates of these models to compute the direct (Figure 1: path c) and indirect effects (joint estimate for Figure 1: paths a and b) between Neuroepo dosage and latent cognitive variable. Estimates for the direct and indirect path are calculated via the mediation package (Tingley et al., 2019) in R with 1000 Monte-Carlo simulations.

Table 1 shows the demographic and clinical characteristics of the sample.

The average age of patients in the Neuroepo group was 56.4 years old (SD = 7.8), and the average age for the placebo group was 60.9 years old (SD = 7.2). The average duration of the disease for patients in the two groups was 5.4 years (SD = 3.2) and 5.8 years (SD = 4.1), respectively. Using the Hoehn and Yahr scale, we classified four patients in the Neuroepo group and one in the placebo group at stage I. We also identified eleven patients from the Neuroepo group and ten in the placebo group at stage II. Only nine patients had familial antecedents of PD (six patients in the Neuroepo group and three in placebo). No significant differences were found between the two groups concerning age (p = 0.16), duration of illness (p = 0.98), familial antecedents of PD (p = 0.5), or severity of PD (p = 0.26.). Additionally, we found no significant differences between baseline and 6-months after intervention in both groups for the total score of the UPDRS “on” (p = 0.71) and “off” (p = 0.88), neither the LED (p = 0.77). There were no significant differences in mood between both groups at baseline.

The whole qEEG spectra of (3,244 voxels × 49 frequency points) for each subject were used to compute the pCCA whitening transform. For this analysis, the six cognitive scores selected: (DRS, FAB, MMSE, Sequency, Memory, and Recognition) were included in calculating the latent cognitive scores. pCCA whitening was carried out to find a joint transform for the cognitive scores and the complete set of qEEG variables.

The pCCA of these two datasets showed that the first latent dimension shows the strongest positive association between qEEG and Cognition. We restricted the further analysis to the first dimension, which was most highly associated with a canonical correlation of 0.97.

Table 3 and Figures 2, 3 show the loadings of the first dimension from CCA whitened qEEG and Cognition.

The Cognitive latent variable is shown in Table 3 shows the loading of W_Cog. The first column is the abbreviation of the cognitive variable and its respective loadings in the second column. Note that the loadings are roughly the same value. Thus, this latent variable is approximately an average of all the tests. The highest loadings are for DRS, followed by FAB, while the minimum loading is for Memory within this narrow range.

The marginal distributions of the loadings for the qEEG latent variables W_qEEG are shown in Figure 2 (frequency) and Figure 3 (space). Figure 2 shows the frequency spread of the latent variable by displaying the minimum and maximum loadings for the 3,244 sources at each of the frequency points. In Figure 2, the x-axis is the frequency in Hz with dashed vertical lines for each classical EEG band (delta, theta, low alpha, high alpha, and beta), and the y-axis is the standardized latent loadings. The red line shows the highest loadings, and the cyan line shows the minimum loadings for each frequency point. The latent factor maxima and minima are loads distributed over all frequency bands. This pattern suggests that this latent variable is a contrast (weighted linear combination) between different cortical areas at each frequency.

This impression is substantiated by Figure 3 shows the dominant loadings (top 1% positive and negative values) plotted on the cortex and summarized for the classical frequency bands. There is a complex pattern of positive and negative values.

We have implemented two mixed effect models with the latent variables for cognition and qEEG (mediator and outcome models).

The mediator model showed that qEEG depends positively on dose and education with a p-value of 0.003 and 0.02, respectively. In contrast, the fixed effects for age showed a negative coefficient; however, the p-value was not significant. Results are given in Table 4. The mediation model’s power analysis showed that the dosage effect is sufficiently powered with 88.00% (95% CI [79.98, 93.64]).

The second model (outcome model) for Cog shows a positive dependency on qEEG and dose with a p-value of <0.0001 and 0.006, respectively (Table 5). The effect of age on Cognition also showed a negative trend; however, the p-value was not significant, possibly due to the small sample size and narrow age range. We performed the power analysis for the outcome model as well. The qEEG and dose estimate had sufficient power of 100.0% (95% CI [96.38, 100.0]) and 89.00% (95% CI [81.17, 94.38]).

Moreover, we combined these two models using the mediation package. We computed the estimates and quasi-Bayesian confidence intervals for the mediation models (Figure 1), the direct effect pathway (Figure 1: path c), and the mediation effect pathways (Figure 1: paths a and b).

The results showed a positive association between dose and Cognition, meaning a higher value for dose results in higher values for the latent cognitive scores (Cog). The beta estimate for the total effect was 2.9 with 95% CI [1.7–4.1] and p-value < 0.0001. The 66% of the total effect was mediated by qEEGwith a p-value < 0.0001. The mediation analysis results are summarized in Table 6 and Figure 4. There is also a significant direct effect between dose and Cog (p-value 0.002), which is not mediated by qEEG. Figure 4 shows the point estimates and 95% CI for mediation (indirect), direct, and total effect.

First, the clinical trial was designed to test Neuroepo tolerance, where safety and not efficacy were the primary outcome measure. However, phase I studies can have therapeutic potential, mainly in patients and not healthy volunteers (Horton, 2006).

Second, the generalizability of our findings is limited because of the size of the sample. The phase I–II trials use small samples since they are testing novel drugs. However, the output was sufficiently powered. For example, qEEG, when analyzed with dose, reached 89% (79.9, 93.6), and analyzed with Cognition reached 100% (96.3, 100).

Third, with such a short treatment and follow-up, we cannot entirely rule out that an underlying idiopathic change in PD progression causes the observed effects on Cognition. The LME contrast would eliminate a factor common to all subjects. This process would have to be distributed differentially in both groups. Larger sample sizes, longer follow-ups, and additional biological measurements would be needed to identify this process.