Verónica Muñoz-Ramírez

Nicolas Pinon

Carole Lartizen

Michel Dojat

Patch vs. global image-based unsupervised anomaly detection in MR brain scans of early Parkinsonian patients

Keywords: Parkinson's disease, Anomaly detection, Patches, Siamese networks, Auto-encoders

Although neural networks have proven very successful in a number of medical image analysis applications, their use remains difficult when targeting subtle tasks such as the identification of barely visible brain lesions, especially given the lack of annotated datasets. Good candidate approaches are patch-based unsupervised pipelines which have both the advantage to increase the number of input data and to capture local and fine anomaly patterns distributed in the image, while potential inconveniences are the loss of global structural information. We illustrate this trade-off on Parkinson's disease (PD) anomaly detection comparing the performance of two anomaly detection models based on a spatial auto-encoder (AE) and an adaptation of a patch-fed siamese auto-encoder (SAE). On average, the SAE model performs better, showing that patches may indeed be advantageous.

Introduction

Medical imaging represents the largest percentage of data produced in healthcare and thus a particular interest has emerged in deep learning (DL) methods to create support tools for radiologists to analyze multimodal medical images, segment lesions and detect subtle pathological changes that even an expert eye can miss. The vast majority of these methods are based on supervised models which require to be trained on large series of annotated data, time-and resourceconsuming to generate.

Over the years, several publicly available neuroimaging databases have been curated and completed with annotations. Some of the most prominent ones are: MSSEG, for multiple sclerosis lesion segmentation [START_REF] Commowick | Objective evaluation of multiple sclerosis lesion segmentation using a data management and processing infrastructure[END_REF]; BRATS, for brain tumor segmentation [14]; ISLES, for ischemic stroke lesion segmentation [START_REF]MICCAI: Ischemic stoke lesion segmentation challenge[END_REF]; and mTOP for mild traumatic injury outcome prediction [START_REF]MICCAI: Mild traumatic brain injury outcome prediction[END_REF]. Challenges, namely those of MICCAI, are organized regularly to showcase the latest technological advancements and push the community towards better performances. However, there are several neurological diseases seldom studied due to the small size and subtlety of the lesions they present. This is the case of vascular disease, epilepsy, and most neurodegenerative diseases in their early stages. The main challenge for such pathologies indeed, is to identify the variability of the pathological patterns on images where the lesion is barely seen or not visible.

Unsupervised methods are good candidates to tackle both the lack of annotated examples and the subtlety of brain scan anomalies [START_REF] Baur | Autoencoders for unsupervised anomaly segmentation in brain MR images: A comparative study[END_REF][START_REF] Shinde | Predictive markers for Parkinson's disease using deep neural nets on neuromelanin sensitive MRI[END_REF]. They rely on networks that learn to encode normal brain patterns in such a manner that any atypical occurrence can be identified by the inability of the network to reproduce it. Auto-encoders (AE), variational auto-encoders (VAE) [START_REF] Kingma | Auto-Encoding Variational Bayes[END_REF] and generative adversarial networks (GAN) [START_REF] Goodfellow | Generative Adversarial Networks[END_REF], have been extensively used as building blocks for unsupervised anomaly detection due to their ability to learn high-dimensional distributions [START_REF] Baur | Autoencoders for unsupervised anomaly segmentation in brain MR images: A comparative study[END_REF].

Parkinson's Disease (PD) is a neurodegenerative disorder that is only identifiable through routine MR scans at an advanced stage. Nevertheless, the manifestation of non-motor symptoms, years before the apparition of the first motor disturbances, suggests the presence of physio-pathological differences that could allow for earlier PD diagnosis. PD afflicts patients for as many as one to two decades of their lives and current treatments can only attenuate some motor manifestations [START_REF] Zhao | Progression of Parkinson's disease as evaluated by Hoehn and Yahr stage transition times[END_REF]. Therefore, reducing the gap between diagnosis and the onset of the neurodegenerative process is of paramount importance to identify personalized treatments that would significantly slow its natural progression. Unsupervised anomaly detection models are here employed to explore such challenging MR data analysis.

In a previous work [START_REF] Muñoz Ramírez | Deep learning models to study the early stages of parkinson's disease[END_REF], we compared deterministic and variational, spatial and dense autoencoders for the detection of subtle anomalies in the diffusion parametric maps of de novo (i.e., newly diagnosed and without dopaminergic treatment) PD patients from the PPMI database [START_REF] Marek | The parkinson's progression markers initiative (ppmi) -establishing a pd biomarker cohort[END_REF]. Our results, while preliminary, offered compelling evidence that DL models are useful to identify subtle anomalies in early PD, even when trained with a moderate number of images and only two parametric maps as input.

Our goal in this paper is to compare an improved anomaly detection pipeline based on a deterministic spatial auto-encoder, hereafter simply referred as AE, to an adaptation of patch-based siamese auto-encoder (SAE) proposed in [START_REF] Alaverdyan | Regularized siamese neural network for unsupervised outlier detection on brain multiparametric magnetic resonance imaging: Application to epilepsy lesion screening[END_REF]. This architecture was originally intended to the detection of subtle epileptic lesions, application for which it achieved promising results.

One important difference between the two compared architectures is the dimension of the input and output data. While AE were trained on 2D transverse slices, thus capturing a global pattern in the image, SAE were trained on small patches sampled throughout the data, making them more suitable to capture fine patterns but losing global structural integrity. Through this comparison we aim to analyze the advantages of patch-fed architectures for the identification of subtle and local abnormalities as well as to propose an alternative for anomaly detection in moderate size image datasets.

Brain anomaly detection pipeline

The anomaly detection task with auto-encoders can be formally posed as follows:

-An auto-encoder is first trained to reconstruct normal samples as accurately as possible. This network is composed of two parts: an encoder (1) that maps the input data into a lower dimensional latent space, assumed to contain important image features, and a decoder (2) that maps the code from the latent space into an output image. -When fed by an unseen image, this trained network produces a reconstructed image from its sampled latent distribution which is the counterpart 'normal' part of the input image. -Reconstruction error maps, computed as the difference between the input and output images, are thus assumed to highlight anomalous regions of the input data.

-Anomaly scores at the voxel, region of interest or image levels may then be derived from the post-processing of these reconstruction error maps.

In this work, we present a general framework for unsupervised brain anomaly detection based on auto-encoders to produce reconstruction error maps and a novel post-processing step to derive per-region anomaly scores.

Autoencoder architectures

We constructed and evaluated two auto-encoder models: a classic auto-encoder (AE) and a siamese auto-encoder (SAE). Both models are fully-convolutional. Their architectures are displayed in Figure 1 and their differences are detailed below.

Classic Auto-encoder: This architecture consists of 5 convolutional layers that go from input to bottleneck and 5 transposed convolutional layers going from bottleneck to output. The output of the encoder network is directly the latent vector z and the loss function was simply the L 1 -norm reconstruction error between input x and output x: The second term of the loss function L SAE (eq. 2) is designed to maximize the cosine similarity between z 1 and z 2 . This constraints patches that are "similar" to be aligned in the latent space. Unlike standard siamese architectures where similar and dissimilar pairs are presented to the network, Alaverdyan et al. [START_REF] Alaverdyan | Regularized siamese neural network for unsupervised outlier detection on brain multiparametric magnetic resonance imaging: Application to epilepsy lesion screening[END_REF] proposed to train this architecture on similar pairs only and compensate for the lack of dissimilar pairs through a regularizing term that prevents driving the loss function to 0 by mapping all patches to a constant value. This term is defined as the mean squared error between the input patches and their reconstructed outputs. The proposed loss function for a single pair hence is:

L AE (x) = x -x 1 (1)
L SAE (x 1 , x 2 ) = 2 t=1 ||x t -xt || 2 2 -α • cos(z 1 , z 2 ) ( 2
)
where xt is the reconstructed output of the patch x t while z t is its representation in the middle layer bottleneck and α an hyperparameter that controls the tradeoff between the two terms. As depicted on Figure 1, the encoder part is composed of 3 convolutional layers and one maxpooling layer in-between the first and second convolutions, while the (non symmetrical) decoder part is composed of 4 convolutional layers, with an upsampling layer in-between the second and third convolutional layer.

Post-processing of the reconstruction error maps

We leveraged the reconstruction error maps obtained from both architectures to generate an anomaly score, following the methodology introduced in [START_REF] Muñoz-Ramírez | Subtle anomaly detection in mri brain scans: Application to biomarkers extraction in patients with de novo parkinson's disease[END_REF]. The voxel-wise reconstruction errors in one image were computed as ||x i -xi || 1 . Since the architectures were fed more than one channel (here two MR modalities), we defined the joint reconstruction error of every voxel as the square root of the sum of squares of the difference between input and output for every channel.

Next, we fixed a threshold on these generated reconstruction error maps to decide whether or not a given voxel should be considered as abnormal, hereafter called the abnormality threshold. Since we expected PD patients to exhibit abnormal voxels in larger quantities than controls, this value corresponded to an extreme quantile (eg. the 98% quantile) of the errors distribution in the control population. The thresholded reconstruction error maps were then employed to identify anatomical brain regions for which the number of abnormal voxels could be used to discriminate between patients and controls.

Experiments

Data

The dataset used in this work consisted of DTI MR scans of 57 healthy controls and 129 de novo PD patients selected from the PPMI database. All images were acquired with the same MR scanner model (3T Siemens Trio Tim) and configured with the same acquisition parameters. Only one healthy control was taken out of the study due to important artifacts in the images.

From these images, two measures, mean diffusivity (MD) and fractional anisotropy (FA), were computed using MRtrix3.0. Values of FA and MD were normalized into the range [0, 1]. The images were spatially normalized to the standard brain template of the Montreal Neurological Institute (MNI) with a non-linear deformation. The resulting MD and FA parameter maps were of dimension 121 × 145 × 121 with a voxel size of 1.5 × 1.5 × 1.5mm 3 .

The control dataset was divided into 41 training controls and 15 testing controls to avoid data leakage. This division was effectuated in 10 different manners through a bootstrap procedure in order to assess the generalization of our predictions as advised in [START_REF] Poldrack | Establishment of Best Practices for Evidence for Prediction: A Review[END_REF]. We took special care to maintain an age average around 61 years old for all the training and test population as well as a 40-60 proportion of females and males.

Once the models were trained with one of the 10 training sets, they were evaluated with the corresponding healthy control test set and the PD dataset (age: 62 y. ± 9; sex: 48 F).

Training of the auto-encoders AE models

The training dataset of the AE models consisted of 1640 images corresponding to 40 axial slices around the center of the brain for each of the 41 training control subjects. The AE models were trained for 160 epochs, with a learning rate of 10 -3 . 3 × 3 kernels were convoluted using padding of 1 pixel and a stride of (2, 2). The bottleneck dimensions were h=4, w=5 and c=256. There were no pooling layers. Implementation was done in Python 3.6.8, Py-Torch 1.0.1, CUDA 10.0.130 and trained on a NVIDIA GeForce RTX 2080 Ti GPU with batches of 40 images. After each convolutional layer, batch normalization [START_REF] Ioffe | Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[END_REF] was applied for its regularization properties. The rectified linear unit (ReLU) activation function was employed in each layer except the last, for which a sigmoïd was preferred. The loss functions were optimized using Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF].

SAE models

The SAE model was trained with 600 000 patches of size 15×15×2 (∼15 000 patches per subject). The model was trained for 30 epochs, with a learning rate of 1 × 10 -3 . Bottleneck dimensions were h=2, w=2 and c=16. Maxpooling and upsampling layers were used, as detailled before. No batch-normalization layers were used. Activation function for every convolutional layer was the rectified linear unit (ReLU) and the sigmoid function was used in the last layer. The kernel size and the numbers of filters were 3x3 and 16 respectively for all convolution blocks but the final one with 2x2 and 2 (equal to the number of channels) respectively. The stride was 1 for all convolution blocks. The maxpooling/upsampling factor was 2. Implementation was done in Python 3.8.10, Tensorflow 2.4.1, 11.0.221 and trained on a NVIDIA GeForce GTX 1660 GPU with batches of 225 patches. The loss function was comprised of a reconstruction part (mean squared error) and a similarity measure (cosine similarity) weighted by a coefficient α=0.005. The loss function was minimized using Adam [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF].

Performance evaluation

The percentage of abnormal voxels found in the thresholded reconstruction error maps was employed to classify them as healthy or pathological (PD) based on a threshold. The critical choice of the threshold was investigated using a Receiver Operator Curve (ROC), taking into account the imbalanced nature of our test set (15 healthy and 129 PD). Every point in the ROC corresponds to the sensitivity and specificity values obtained by a given threshold. As proposed in [START_REF] Muñoz-Ramírez | Subtle anomaly detection in mri brain scans: Application to biomarkers extraction in patients with de novo parkinson's disease[END_REF], the choice of the optimal threshold, referred to as the pathological threshold, was based on the optimal geometric mean, g-mean= √ Sensitivity × Specif icity. Additionally, to help evaluate the localization of anomalies, two atlases were considered: the Neuromorphometrics atlas [START_REF] Bakker | The Scalable Brain Atlas: Instant Web-Based Access to Public Brain Atlases and Related Content[END_REF] and the MNI PD25 atlas [START_REF] Xiao | Multi-contrast unbiased MRI atlas of a Parkinson's disease population[END_REF]. The first was used to segment the brain into 8 macro-regions: subcortical structures, white matter and the 5 gray matter lobes (Frontal, Temporal, Parietal, Occipital, Cingulate/Insular). The latter was specifically designed for PD patients exploration. It contains 8 regions: substantia nigra (SN), red nucleus (RN), subthalamic nucleus (STN), globus pallidus interna and externa (GPi, GPe), thalamus, putamen and caudate nucleus. For all of the before-mentioned regions of interest (ROI), we calculated the g-mean for the associated pathological threshold, leading to the classification performance of our models.

Results

As it can be seen in Figure 2, both auto-encoder architectures achieve good quality reconstructions, however the SAE seems to capture finer details and textures than the AE. This explains the high contrast in AE reconstructions. The visualization of the percentage of abnormal voxels in the ROIs presented in Figure 3 showcases the inter-subject variability amongst members of the same population (healthy and PD). Even so, abnormal voxels are clearly more numerous in the PD patient population. The g-mean classification scores for all models, obtained for each ROI and each sub-population, are presented in Figure 4. We notice that, on average, the SAE model performed better than the AE on the whole brain and most of the macro and subcortical structures studied, with the exception of the temporal lobe, the putamen, the thalamus and the internal and external segment of the globus pallidus. We note that the results varied greatly across the ten populations samples. As an indication, for the whole brain, the SAE obtained a g-mean average score of 66.9 ± 5.8% and the AE 65.3 ± 7.5%, however the best scores among the 10 samples were of 79.9% and 81.9% for the AE and SAE respectively, both on sample 1. Corresponding values for the white matter are of 68.2 ± 4.6% for SAE and 66.2±6.7% for AE. The largest standard deviations in the observed anatomical regions belonged to the white matter, the frontal and occipital lobes.

Discussion and conclusion

Unsupervised auto-encoders (AE) have shown to efficiently tackle challenging detection tasks where brain alteration are barely seen or not visible. The objective of our study was to explore the potential of such AE models for the detection of subtle anomalies in de novo PD data and compare patch-based versus image-based models.

Both AE and SAE architectures produced good quality reconstructions and were able to discriminate between healthy individuals and recently diagnosed PD patients with performances (see Figures 3 and4) that are competitive with those found in the literature. Notably, the Correia & al. [START_REF] Correia | Towards accurate and unbiased imaging-based differentiation of Parkinson's disease, progressive supranuclear palsy and corticobasal syndrome[END_REF] SVM mean accuracy score for a selection of WM regions is of 61.3% whereas both the SAE and AE achieved g-mean scores above 66% for the WM. Also, the cross-validation procedure of Schuff & al. [START_REF] Schuff | Diffusion imaging of nigral alterations in early Parkinson's disease with dopaminergic deficits[END_REF] obtained a ROC AUC of 59% for the rostral segment of the SN which is below our SN average g-mean score for the SAE and equal to that of the AE.

Note that at this early stage of the disease (1-2 H-Y scale) the patients have no tremor nor uncontrolled movements compared to healthy subjects. This rules out that movement was the index that allowed PD classification.

Using DTI data we did not search for structural atrophy or lesion load but for degradation of WM properties in the early stages of PD that could appear everywhere in the brain. This explains why the WM obtains the highest gmean scores. This being said, our models could largely improve by increasing the size of our dataset. Furthermore, the addition of another MR modality such as iron load using T2/T2* relaxometry could allow us to detect the reduction of dopaminergic neurons in subcortical structures, largely reported in the early stages of PD but not visible in DTI.

Regarding the comparison between the AE and SAE, the choice is not clear. While the classic AE architecture benefits from a more straightforward implementation, the SAE proposes significant advantages for small databases. Indeed, patch-fed networks can be trained with smaller samples of data and the siamese constraint of the architecture ensures efficient learning. What is more, the latent space features of these models contain local information that can be used to classify between healthy and pathological individuals at the voxel-level and produce anomaly maps like those introduced by [START_REF] Alaverdyan | Regularized siamese neural network for unsupervised outlier detection on brain multiparametric magnetic resonance imaging: Application to epilepsy lesion screening[END_REF].

In future work, we plan to generate this kind of maps for early PD patients to offer more precise indications about the localization of anomalies and correlate them with the PD hemispheric lateralization. In addition, we aim to complete our dataset by adding other MR modalities such as perfusion and relaxometry, but also by gathering heterogeneous data from multi-vendors and multi-sites exams. For this purpose, we will use an harmonization procedure as a preprocessing step (extension of DeepHarmony [START_REF] Dewey | DeepHarmony: A deep learning approach to contrast harmonization across scanner changes[END_REF]). Finally, our 3D implementation of the SAE model is ongoing.

Fig. 1 .

 1 Fig. 1. Classical auto-encoder (AE) on top, Siamese auto-encoder (SAE) at the bottom

Fig. 2 .

 2 Fig. 2. Showcase of a slice of the original data and its AE and SAE reconstructions

Fig. 3 .

 3 Fig. 3. The percentage of abnormal voxels found by the SAE in the anatomical ROIs presented in Section 3. Top: the test controls of Sample 1; Bottom: 15 randomly selected PD patients.

Fig. 4 .

 4 Fig. 4. g-mean scores for the whole brain and several ROIs for AE and SAE. The vertical dashed lines separate macro-and micro brain structures.


Duje Buric

Christine Chacon

Gerard Audran

Diane Braguer

R A Sylvain

Manon Marque

Paul Carre

Toshihide Bremond

Yamasaki

Toshihide Yamasaki

Gérard Audran

Marque

Manon Carré
email: manon.carre@univ-amu.fr

Paul Brémond
email: paul.bremond@univ-amu.fr

Chemical modifications of imidazole-containing alkoxyamines increase C-ON bond homolysis rate: Effects on their cytotoxic properties in glioblastoma cells

Keywords: Alkoxyamines, Alkyl radicals, Glioblastoma, Antitumor activity

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The blood-brain barrier (BBB) is a semipermeable membrane barrier indispensable for the normal function of the central nervous system (CNS). Not only it hinders influx of most compounds and microbes from blood to brain but also it blocks most of therapeutics from reaching their targets in the brain. As a consequence, very few chemotherapeutic agents are available for treatment of brain tumors, such as glioblastoma, which are among the most aggressive and common brain tumor in adults as well as children. In the absence of active transport through the BBB, the molecular weight (MW), lipophilicity (logP), distribution coefficient (logD), hydrogen bond donors (HBD) and ionization of compound (pKa) are among the most important factors for the diffusion of the therapeutic agents [1]. Therefore, novel low-MW molecules active against cancer cells are exceptionally interesting in the field of cancer therapeutics since they are more prone to penetrate BBB and thus might be used against glioblastoma.

Alkoxyamines, with R 1 -ONR 2 R 3 general structure, are compounds able to undergo homolysis of the labile C-ON bond to release a transient alkyl radical (R 1 •) and a generally stable nitroxide (R 2 R 3 NO•) [2][3][START_REF] Prescott | Biological Relevance of Free Radicals and Nitroxides[END_REF]. The main application of these compounds is the nitroxidemediated controlled radical polymerization [5]. Interestingly, we recently proposed a new application of these low-molecular weight compounds as theranostic agents, i.e. with both therapeutic and diagnostic roles in fighting malignant tumors [6]. Indeed, we demonstrated that the released nitroxide from a pyridine-containing alkoxyamine can be used to perform imaging by the Overhauser-enhanced MRI (OMRI) technique, while the released alkyl radical triggered the apoptosis of cancer cells in vitro (Scheme 1) [START_REF] Moncelet | Alkoxyamines: Toward a New Family of Theranostic Agents against Cancer[END_REF].

Scheme 1 Alkoxyamines as theranostic agents

However, the low pKa of the pyridine ring of alkoxyamine 1a (pKa = 4.70) was an unfavorable aspect since it precludes the protonation of the ring in physiological conditions, i.e. at pH 7.4, and it required a beforehand activation by methylation to obtain the labile alkoxyamine 1b (t1/2 = 45 min) [8]. Furthermore, 1b showed moderate cytotoxicity against glioblastoma cell lines. Therefore, these results prompted us to replace the pyridine ring of 1a by an imidazole ring (pKa = 6.95), a moiety that is widely considered as biocompatible and is present in many natural compounds and drugs [9]. Herein, we report the synthesis of sixteen new alkoxyamines that contain an imidazole ring grafted at the C-2 (2a-c), N-1 (3a-b) or C-4 (4a-c) positions, respectively, and in this study, we will determine their kinetic properties and assess their therapeutic potential for glioblastoma (Figure 1). All these alkoxyamines are based on the same nitroxide part as 1b, i.e. SG1 nitroxide, as it affords labile alkoxyamines and as its imaging properties were already highlighted [START_REF] Moncelet | Alkoxyamines: Toward a New Family of Theranostic Agents against Cancer[END_REF].

Results and discussion

Synthesis of new alkoxyamines containing an imidazole ring

Synthesis of an alkoxyamine generally requires the in-situ generation of an alkyl radical followed by nitroxide-trapping of the latter. One of the most convenient methods for this purpose is the combination of an activated double bond and a particular catalytic system, based on a manganese-salen complex, with NaBH4 in i-PrOH and under air atmosphere. Therefore, Mn-catalyzed coupling of SG1 nitroxide (N-(2-methylpropyl)-N-(1diethylphosphono-2,2-dimethylpropyl)-N-oxyl) and 2-vinyl-1H-imidazole gave alkoxyamine 2a in 54% yield and as two diastereoisomers, which were purified and separated by column chromatography to give (RR/SS)-2a and (RS/SR)-2a (Scheme 2). The same procedure was repeated to obtain alkoxyamines 3a (37% yield) and 4a (39% yield) from 1-vinyl-1Himidazole and 4-vinyl-1H-imidazole, respectively, as well as mono-methylated analogue 2b (54% yield) from 1-methyl-2-vinyl-1H-imidazole. Dimethylimidazolium analogue 2c and monomethylimidazolium analogue 3b were obtained by methylation of each diastereoisomer

N O P(O)(OEt) 2 N O P(O)(OEt) 2 N N N N N N N N N N N N (RR/SS)- and (RS/SR)- diastereoisomers 2a 2b 2c 4a 4b 4c H H N N N N 3a 3b
methylated analogue 4b and dimethylimidazolium analogue 4c were both obtained directly by methylation of 4a using methyl iodide and potassium carbonate in dichloromethane (see experimental section). All compounds were obtained as pure and diastereomerically homogeneous compounds and were characterized by NMR ( 1 H, 13

Kinetic investigations

The pKa values of alkoxyamines were determined by titration by 1 H NMR of signals of the aromatic moiety at various pH in D2O/CD3OD (1:1 v/v) [11]. All pKa were ca. 1-2 units higher than the previously investigated alkoxyamine 1a, that contains a pyridine ring (Table 1).

Then, the stability of the new alkoxyamines was investigated by EPR in aqueous solutions, using oxygen as radical scavenger. Thermolysis experiments were performed at different pH,

i.e. to obtain >99% of protonated or non-protonated forms of imidazole rings. Homolysis rates (kd) were obtained and the corresponding activation energies (Ea) were calculated using Arrhenius equation and the averaged pre-exponential factor A = 2.4.10 14 s -1 . Half-life times (t1/2) were also re-estimated at 37 ºC to allow direct comparisons between compounds at physiological temperature. Depending on pH, clear differences in Ea and t1/2 values were observed (Table 1 and supplementary material). values similar to other SG1-based alkoxyamines bearing a secondary benzylic-type alkyl radical. Small differences of 0.3-0.7 kJ/mol were observed between diastereoisomers within the conventional ranges, as already documented [START_REF] Audran | C-ON bond homolysis in alkoxyamines. Part 12: the effect of the para-substituent in the 1-phenylethyl fragment[END_REF]. From these Ea, t1/2 of 10-19 d at 37 ºC were calculated. However, (RR/SS)-and (RS/SR)-3a were found to be very stable in aqueous solutions. Therefore, DMSO was used as solvent in order to reach higher temperature and accelerate thermolysis experiments. Very high Ea were determined and, consequently, estimated t1/2 at 37 ºC clearly illustrated the high stability of alkoxyamine 3a. The differences of stability between the alkoxyamines were explained by performing calculations of the radical stabilization energy (RSE) of benzylic radical models. Indeed, the RSE of the calculated radicals were -3.6 kJ/mol for the series 2 and +4.3 kJ/mol for the series 4 (see supplementary material). These values indicated that the stabilization of the odd electron is better in the series 2 than in the series 4, thus confirming the observed reactivity. On the other hand, in the case of the series 3, the calculated RSE was highly positive (+20.4 kJ/mol), i.e.

small or no delocalization of the odd electron, thus confirming the exceptional stability for this SG1-based alkoxyamine [START_REF]Ea values of (RR/SS)-and (RS/SR)-diastereoisomers of the benchmark secondary alkoxyamine SG1-CH(Ph)CH3 are 123.0 and 124.0 kJ/mol, respectively[END_REF].

At acidic pH, as expected, acceleration of the homolysis rate by protonation of the imidazole ring was observed for (RR/SS)-2a and (RS/SR)-2a, as well as (RR/SS)-4a and (RS/SR)-4a.

However, the efficiency of the proton-triggered activation strongly depends on the position of the nitrogen atoms in the imidazole ring. Indeed, we observed decreases in Ea of ca. 10.5 kJ/mol for (RR/SS)-2a and (RS/SR)-2a, while it was only ca. 2-4 kJ/mol for (RR/SS)-4a and

(RS/SR)-4a. For 3a, upon protonation with PTSA in DMSO, activation was also observed, as

Ea of (RR/SS)-3a decreased of 6.1 kJ/mol, while it decreased of 4.0 kJ/mol for (RS/SR)-3a.

Such differences in reactivity depending on the position of the nitrogen atoms were already observed by us on pyridine-containing alkoxyamines 1a and other analogues [START_REF] Audran | Chemically triggered C-ON bond homolysis in alkoxyamines. Part 7. Remote polar effect[END_REF]. Indeed, while the development of the positive charge due to protonation affects the entire heteroaromatic ring, thus increasing its polar effect, it is however more prevailing on particular positions of the imidazolium ring, i.e. N-1, C-2 and N-3 positions (Scheme 3).

Hence, in the series 2, a positive charge can be developed on the aromatic carbon atom in αposition of the labile bond (i.e. C-2 position), thus strong chemical activation of the homolysis, i.e. acceleration of 60-times, can be achieved. On the other hand, in the series 4, the positive charge cannot be developed in the adequate position to obtain chemical activation, which would be the C-4 position, thus only the increase of the polarity, due to the protonation of the heteroaromatic moiety, is probably the sole responsible for the activation, with only a small effect (ca. 2-4-times acceleration). In the series 3, despite the apparent favorable development of the positive charge in the adequate position (i.e. N-1 position), the strong destabilization of the released radical leads to a high stability of the corresponding alkoxyamines (vide supra) [START_REF] Moreover | anomeric effect, i.e. nO → σ*C-N, that may significantly strengthen the C-ON bond, was observed on the X-ray structure of (RS/SR)-3a, nor on the calculated structures[END_REF].

Scheme 3 Mesomeric forms of alkoxyamines (R, R' = H or Me)

N N R R' SG1 N N R R' SG1 N N R R' SG1 N N SG1 N N SG1 N N SG1 R' R R R' R R' N N SG1 N N SG1 N N SG1 R R R 3a-b 2a-c

4-ac

Monomethylation of 2a into 2b afforded a loss of 1.6 to 5.4 kJ/mol at basic pH depending on the diastereoisomers, i.e. acceleration of homolysis of 2-to 8-times. After protonation, an additional decrease of 7.7 to 10.3 kJ/mol, i.e. acceleration of 20 to 50-times, was observed.

Thus, by combining methylation and protonation, the total decrease was ca. 12-13 kJ/mol, meaning ca. 100-to 150-times acceleration of homolysis rate. This increased activation can be ascribed to an increase of the steric hindrance of the alkyl part, leading to an acceleration of the homolysis rate. However, methylation of 4a into 4b was less striking, leading to an overall decrease of ca. 3-5 kJ/mol. i.e. only ca. 3-to 6-times acceleration, probably because the increase of the steric hindrance due to the methylation is farther away from the reactive center and therefore its effect is negligible.

In order to grasp their potential reactivity in biological conditions, kinetics of all these alkoxyamines were also measured at pH 7.4 (Table 2). Due to its very high Ea, kinetic of alkoxyamine 3a was not re-measured. However, intermediates values of Ea, as well as of t1/2, were obtained, as pKa values of these compounds are lower than 7.4, thus limiting the activation of the homolysis. Therefore, in order to circumvent this issue and to potentially achieve full activation at pH 7.4, we followed the same approach used for alkoxyamine 1a by converting the imidazole ring into the corresponding imidazolium by methylation investigating dimethylimidazolium-alkoxyamines 2c, 3b and 4c at pH 7.4. For alkoxyamine 2c, overall decreases of 8.1 to 11.7 kJ/mol, i.e. acceleration of 24-to 90-times, were observed, thus confirming that it is unstable at pH 7.4. Similarly, for alkoxyamine 4c, a decrease of Ea (ca. 7 kJ/mol) was observed, i.e. acceleration factor of ca. 15-times. However, for alkoxyamine 3b, despite that decreases of ca. 10 to 14 kJ/mol were observed, half-life times were still very high. b Calculated with experimental kd values using Arrhenius equation and A = 2.4.10 14 s -1 .

c Calculated at 37 ºC using re-estimated kd values.

d Not determined (see text).

Cytotoxic properties of the alkoxyamines

We previously reported that alkoxyamines display lethal properties in glioblastoma cells, but with a poor level of efficacy [START_REF] Moncelet | Alkoxyamines: Toward a New Family of Theranostic Agents against Cancer[END_REF]. Therefore, we investigated whether the cytotoxicity of the new unstable alkoxyamines (2a-c and 4a-c) was improved in human glioblastoma cell lines and whether stable alkoxyamines (3a-b) stay inert due to their extremely long half-life. The cytotoxic potential of the sixteen novel alkoxyamines (2a-c, 3a-b and 4a-c) was first evaluated in the human U87-MG cell line, derived from malignant glioma that we previously used for research of new active compounds in glioblastomas [START_REF] Le Grand | Akt targeting as a strategy to boost 36 chemotherapy efficacy in non-small cell lung cancer through metabolism suppression[END_REF][START_REF] Denicolaï | Proscillaridin A is cytotoxic for glioblastoma cell lines and controls tumor xenograft growth in vivo[END_REF][START_REF] Alves | A proapoptotic peptide conjugated to penetratin selectively inhibits tumor cell growth[END_REF][START_REF] Mercurio | Evidence for new targets and synergistic effect of metronomic celecoxib/fluvastatin combination in pilocytic astrocytoma[END_REF]. All alkoxyamines exerted a dose-dependent inhibition of the cancer cell survival after a 72 h-treatment (Figure 2a andb). IC50 values, i.e. concentrations that half-reduced U87-MG cell survival, ranged from 25 to 380 µM (Figure 2c). Among the twelve unstable alkoxyamines (2a-c and 4a-c) evaluated in the present study, eight were more effective in reducing cancer cell survival than the previously investigated alkoxyamine 1b, whose IC50 value reached 250 µM in similar experimental conditions [START_REF] Moncelet | Alkoxyamines: Toward a New Family of Theranostic Agents against Cancer[END_REF]. It is especially notable with the compound (RS/SR)-4a, which displayed the highest efficacy and was 10 times more active than 1b in the same U87-MG cell model (Figure 2b andc). The promising anticancer property of this novel alkoxyamine was confirmed by light-microscopy, with a significant decrease in cell density as compared with the other molecules (Figure 2d).

Among the four stable alkoxyamines (3a-b), three of them exerted a significant cytotoxicity with IC50 ranging from 101 to 166 µM, suggesting that either they release alkyl radicals by an alternative mechanism or they are toxic without undergoing homolysis. Therefore, in order to explore whether the novel alkoxyamines -both stable and unstable -exert their cytotoxicity through the release of alkyl radical, we co-treated U87-MG cells with IC50 concentration of both alkoxyamines and the alkyl radical scavenger Troxerutin [6]. The experiment showed that cytotoxicity of all the twelve unstable alkoxyamines (2a-c and 4a-c) was significantly inhibited by Troxerutin confirming that their efficacy comes from the release of alkyl radical (Figure 2d). However, cytotoxicity of the four stable alkoxyamines (3a-b) remained intact despite the presence of the alkyl radical scavenger (Figure 2d). Their mechanism of action is thus independent of homolysis and will undergo further analysis.

Furthermore, cytotoxicity of the twelve unstable alkoxyamines (2a-c and 4a-c) emerges faster than the chemical half-life should have suggested. Chemical half-life of the alkoxyamines in our study was measured in vitro in aqueous solution, which provides anisotropic conditions that are not present in heterogeneous cell culture medium or within the cytoplasm.

Alkoxyamines are likely to obtain other conformations that increase homolysis rate due to biological environment: e.g. various proteins alkoxyamines may be complexed to, biological membranes they may be inserted in, or viscous intracellular cytosol that may result in conformational changes. Details of the mechanism of activation in cellulo are under investigation and will be addressed in the follow-up of the project.

A common feature of the unstable alkoxyamines (2a-c and 4a-c) is that the charged molecules were less cytotoxic towards cancer cells than the neutral ones. To mention but a few examples, the monomethylimidazole alkoxyamines (RR/SS)-4b and (RS/SR)-4b were 4.4 and 4.9 times more active than the dimethylimidazolium alkoxyamines (RR/SS)-4c and (RS/SR)-4c, respectively. Indeed, it is generally accepted that neutral compounds are more readily able to cross non-polar lipid bilayers, unlike charged compounds where this process is energetically disfavored [START_REF] Manallack | The Significance of Acid/Base Properties in Drug Discovery[END_REF]. Thus, despite the fact that the dimethylimidazolium alkoxyamines have lower t1/2, i.e. undergo faster homolysis, the charged heteroaromatic ring may disturb the targeting of the molecules to cell plasma membrane and/or to intracellular compartments. Molecular weight (MW) and distribution coefficient (logD) value have been found to be the most discriminating factors in determining the permeability of drug candidates. For cases where a drug reaches its target locations by passive mechanisms, the ideal logD has been defined to be up to 4 [START_REF] Gleeson | Generation of a set of simple, interpretable ADMET rules of thumb[END_REF], and either >1.3 for compounds with molecular weight below 414 g.mol -1 or >2.4 for higher molecular weights [START_REF] Waring | Defining optimum lipophilicity and molecular weight ranges for drug candidates-Molecular weight dependent lower logD limits based on permeability[END_REF]. We therefore turned our attention to the noctanol/PBS partition coefficient of the different alkoxyamines at physiological pH. As shown in Table 3, the logD7.4 values measured for alkoxyamines 2c, 3b and 4c confirmed that these molecules were not in the optimum range for biological activity. On the other hand, the molecular weights and logD7.4 values of alkoxyamines 2a-b, 3a and 4a-b make them more likely to enter the cells, in agreement with the U87-MG cell survival data. Furthermore, the results obtained in U87-MG cells were additionally validated in the human malignant astrocytoma U251-MG cells. We confirmed that all sixteen alkoxyamines reduced survival of glioblastoma cells in a dose-dependent manner and that the neutral molecules were more active than the charged ones (Figure 3a andb). Importantly, U251-MG cells do not respond to Temozolomide (IC50 TMZ = 455 µM; data not shown), the standard chemotherapy treatment for glioblastoma, pointing out the fact that alkoxyamines may circumvent drug resistance in these refractory cancer cells. Compared with U87-MG cells, U251-MG cells were even more sensitive to alkoxyamines, as shown by the increase in (RR/SS)-2a and

(RS/SR)-2a cytotoxicity by 38% and 36%, respectively (Figure 3c versus Figure 2c). The same applies to (RS/SR)-4a, which was again the most active alkoxyamine evaluated, with an IC50 value reduced to 20 µM (Figure 3c). Higher efficacy in U251-MG cells might be linked to their glycolytic phenotype that, by releasing a greater amount of lactic acid, may result in a faster acidification of the cell microenvironment [START_REF] Li | Differences in Protein Expression between the U251 and U87 Cell Lines[END_REF]. In support to this, we measured the extracellular pH and showed that it remained at 7.5 in the medium sampled from the 48 hcultures of U87-MG cells while it was reduced to 6.8 in U251-MG cell cultures (data not shown). This acidic cell microenvironment may then promote the homolysis of the prodrugs at a faster rate (Tables 1 and2), and faster enough to observe increased cytotoxicities.

Stable alkoxyamines (3a-b) were again more cytotoxic than expected and their efficacy remained intact upon the combined treatment with the alkyl scavenger Troxerutin confirming that, also in U251-MG cells, toxicity of the alkoxyamines 3a-b does not originate from homolysis and alkyl radical release (Figure 3d).

A crucial aspect of glioblastoma treatment concerns the ability of drugs to penetrate the BBB.

This obstacle to drugs entering the CNS results from the tight junctions between the epithelial cells of the BBB and from the efflux pumps residing on astrocyte membranes of the BBB that counteract drug entry. In close relation to the parameters described above for cancer cell membranes permeability, permeability in brain capillaries has been related to drug lipophilicity, and hence to partition coefficient [1]. For the design of compounds where CNS penetration is required, logD7.4 values were indeed given an optimal range between 2.0 and 4.0 [START_REF] Wager | Moving beyond rules: the development of a central nervous system multiparameter optimization (CNS MPO) approach to enable alignment of druglike properties[END_REF], which is the case for the neutral alkoxyamines 2a and 4a (Table 3). The ability of these pro-drugs to be brain-penetrating drugs is also supported by the pKa values, from 6.04 to 6.76 (Table 1), since compounds with pKa above 5.5 have been identified as more likely to be CNS-penetrant [START_REF] Broccatelli | Improving the prediction of the brain disposition for orally administered drugs using BDDCS[END_REF]. Lastly, glioblastoma drug candidates should not be substrates for efflux transporters, such as the P-glycoprotein (P-gp). Considering this active efflux, it has been proposed that compounds are likely to be non-substrates if they have accessible H-accepting atoms ≤ 4, molecular weight < 400 and a pKa < 8 [START_REF] Didziapetris | Classification analysis of Pglycoprotein substrate specificity[END_REF]. Thus, it can be assumed that the active efflux of the alkoxyamines 2a-b and 4a-b from the CNS may be minimized.

Conclusions

Our interdisciplinary work led to the synthesis and characterization of novel imidazole ringcontaining alkoxyamines, and to the analysis of their activity against human cancer cells related to release of free radicals. The understanding of the structure-activity relationship for the unstable alkoxyamines (2a-c and 4a-c), by protonation and/or methylation, highlighted that their biological activities do not only correlate with their half-life time of homolysis, but rather also result from a tight balance between pKa, logD7.4 and molecular weight. Among the set of new alkoxyamines, we identified (RS/SR)-4a as the compound with the highest antitumor activity, including towards glioblastoma cells that are resistant to conventional chemotherapy. Furthermore, its favorable logD7. [START_REF] Prescott | Biological Relevance of Free Radicals and Nitroxides[END_REF] and pKa values make it a robust candidate for blood-brain barrier penetrating therapeutics against brain neoplasia. Thus, our results open new opportunities for future works to optimize the alkoxyamines with increased pharmacological activity and likely as few side effects as possible.

Experimental

Syntheses

General methods

All solvents and reagents were used as received. Routine monitoring of reactions was performed using silica gel 60 F254 TLC plates; spots were visualized using UV light and phosphomolybdic acid solution in EtOH, followed by heating. Purifications by column chromatography were performed with silica gel 60 (230-400 mesh). 1 H, 13 C, and 31 P NMR spectra were recorded in CDCl3 or DMSO-d6 on 300 or 400 MHz spectrometers. Chemical shifts (δ) in ppm are reported using residual non-deuterated solvents as internal reference for 1 H and 13 C-NMR spectra, or using an internal sealed capillary filled with 85% H3PO4 for 31 P-NMR spectra.

Diethyl (1-((1-(1H-imidazol-2-yl)ethoxy)(tert-butyl)amino)-2,2dimethylpropyl)phosphonate ((RR/SS)-2a and (RS/SR)-2a).

To a stirred solution of salen ligand (0.07 g, 0.25 mmol, 0.05 eq.) in i-PrOH (7 mL) was added MnCl2 (0.03 g, 0.25 mmol, 0.05 eq.) in an opened flask. After stirring with air bubbling at rt for 10 min, a solution of SG1 nitroxide (1.47 g, 5.0 mmol, 1.0 eq.) and 2-vinyl-1H-imidazole [START_REF] Lee | Triazine compounds as kinase inhibitors[END_REF] (0.57 g, 6 mmol, 1.2 eq.) in i-PrOH (10 mL) was added, followed by solid NaBH4 (0.19 g, 5 mmol, 1.0 eq.) added in small portions below 25 ºC. The mixture was stirred at rt for 3 h, and then sat. NaHCO3 and CH2Cl2 were added and extracted with CH2Cl2. The combined organic phase was washed with brine, dried with MgSO4, and concentrated to afford a crude product as a 1.

Diethyl (1-(tert-butyl(1-(1-methyl-1H-imidazol-2-yl)ethoxy)amino)-2,2dimethylpropyl)phosphonate ((RR/SS)-2b and (RS/SR)-2b).

To a stirred solution of salen ligand (0.27 g, 1.0 mmol, 0.05 eq.) in i-PrOH (15 mL) was added MnCl2 (0.13 g, 1.0 mmol, 0.05 eq.) in an opened flask. After stirring with air bubbling at rt for 10 min, a solution of SG1 nitroxide (5.89 g, 20 mmol, 1.0 eq.) and 1-methyl-2-vinyl-1H-imidazole [START_REF] Waser | Hydrazines and azides via the metalcatalyzed hydrohydrazination and hydroazidation of olefins[END_REF] (2.49 g, 23 mmol, 1.2 eq.) in i-PrOH (25 mL) was added, followed by solid NaBH4 (0.76 g, 20 mmol, 1.0 eq.) added in small portions below 25 ºC. The mixture was stirred at rt for 6 h, and the solvent was removed in vacuo. Then sat. NaHCO3 and CH2Cl2 were added and extracted with CH2Cl2. The combined organic phase was washed with brine, dried with MgSO4, and concentrated to afford a crude product as a 2.

, J = 139.8 Hz), 69.

dimethylpropyl)phosphonate ((RR/SS)-3a and (RS/SR)-3a).

To a stirred solution of salen ligand (0.27 g, 1.0 mmol, 0.05 eq.) in i-PrOH (15 mL) was added MnCl2 (0.13 g, 1.0 mmol, 0.05 eq.) in an opened flask. After stirring with air bubbling at rt for 10 min, a solution of SG1 nitroxide (5.89 g, 20 mmol, 1.0 eq.) and 1-vinyl-1H-imidazole (2.82 g, 30 mmol, 1.5 eq.) in i-PrOH (25 mL) was added, followed by solid NaBH4 (0.76 g, 20 mmol, 1.0 eq.) added in small portions below 25 o C. The mixture was stirred at rt for 2.5 h, and the solvent was removed in vacuo. Then sat. NaHCO3 and CH2Cl2 were added and extracted with CH2Cl2.

The combined organic phase was washed with brine, dried with MgSO4, and concentrated to afford a crude product as a 1.25: 1 mixture of diastereoisomers ( 31 P-NMR ratio). The crude was separated by flash column chromatography (CH2Cl2: MeOH = 97: 3) to afford (RR/SS)- To a stirred solution of salen ligand (0.23 g, 0.9 mmol, 0.05 eq.) in i-PrOH (15 mL) was added MnCl2 (0.11 g, 0.9 mmol, 0.05 eq.) in an opened flask. After stirring with air bubbling at rt for 10 min, a solution of SG1 nitroxide (5.09 g, 17 mmol, 1.0 eq.) and 4-vinyl-1H-imidazole [START_REF] Cotterill | Thermal 1,3-trityl migrations in Diels-Alder domino reactions of 1-trityl-4-vinyl-1H-imidazoles[END_REF] (1.87 g, 20 mmol, 1.2 eq.) in i-PrOH (25 mL) was added, followed by solid NaBH4 (0.75 g, 20 mmol, 1.2 eq.) added in small portions below 25 o C. The mixture was stirred at rt for 5.5 h, and the solvent was removed in vacuo. Then sat. NaHCO3 and CH2Cl2 were added and extracted with CH2Cl2. The combined organic phase was washed with brine, dried with MgSO4, and concentrated to afford a crude product as a 1: 1 mixture of diastereoisomers ( 31 P-NMR ratio).

3a (1.
The crude was separated by flash column chromatography (CH2Cl2: MeOH = 100: 0 to 9: 1)

to afford (RR/SS)-4a (1. [START_REF] Cotterill | Thermal 1,3-trityl migrations in Diels-Alder domino reactions of 1-trityl-4-vinyl-1H-imidazoles[END_REF]

Figure 1

 1 Figure 1 Structures of the alkoxyamines synthesized and evaluated in this work

Scheme 2

 2 Scheme 2 Synthesis of alkoxyamines (RR/SS)-and (RS/SR)-2a-c, 3a-b, 4a-c: a) MnCl2, N,N′-bis(salicylidene)ethylenediamine (salen), NaBH4,

Figure 2

 2 Figure 2 Cytotoxicity of alkoxyamines 2a-c, 3a-b and 4a-c towards U87-MG human

Figure 3

 3 Figure 3 Cytotoxicity of alkoxyamines 2a-c, 3a-b and 4a-c towards U251-MG human

  dimethyl-1H-imidazol-3-ium 4-methylbenzenesulfonate ((RR/SS)-2c).To a stirred solution of (RR/SS)-2b (0.45 g, 1.1 mmol) in THF (3 mL), a solution of methyl ptoluenesulfonate (1.2 mL, 1.0 M in THF) was added. The mixture was stirred at rt for 1 d under argon, and then the solvent was removed in vacuo. This residue was triturated with pentane/ Et2O (1:1) and the solvent was removed and dried in vacuo to afford (RR/SS)-2c (0.61 g, 92%). pale yellow oil; 1 H-NMR (CDCl3, 400 MHz) δ 0.88 (br. s, 9 H), 1.18 (s, 9 H), 1.33 (t, J = 7.1 Hz, 6 H), 1.76 (d, J = 7.3 Hz, 3 H), 2.33 (s, 3 H), 3.31 (d, J = 25.7 Hz, 1 H), 4.01 (s, 6 H), 3.90 -4.14 (m, 3 H), 4.15 -4.27 (m, 1 H), 5.36 -5.53 (m, 1 H), 7.13 (d, J = 7.8 Hz, 2 H), 7.75 (d, J = 8.2 Hz, 2 H), 7.77 (s, 2 H); 13 C{ 1 H}-NMR (CDCl3, 75 MHz) δ 15.7 (d, J = 6.6 Hz), 16.1 (d, J = 5.5 Hz), 18.5 (s), 20.7 (s), 27.5 (s), 29.4 (d, J = 5.5 Hz), 35.1 (d, J = 4.4 Hz), 35.4 (s), 59.9 (d, J = 7.7 Hz), 61.2 (d, J = 6.6 Hz), 61.8 (s), 69.5 (d, J = 142.0 Hz), 74.6 (br. s), 123.6 (s), 125.4 (s), 128.1 (s), 138.8 (s), 143.0 (s), 144.9 (s); 31 P{ 1 H}-NMR (CDCl3, 162 MHz): δ 23.81; HRMS (ESI-TOF) m/z: [M-TsO -] + Calcd for C20H41N3O4P 418.2829; Found 418.2828. dimethyl-1H-imidazol-3-ium 4-methylbenzenesulfonate ((RS/SR)-2c). (RS/SR)-2c was synthesized in a same procedure as (RR/SS)-2c with (RS/SR)-2b (0.41 g, 1.0 mmol) and methyl p-toluenesulfonate (1.0 mL, 1.0 M in THF) to give (RS/SR)-2c (0.29 g, 48%). white solid; m.p.: 120 o C (dec.); 1 H-NMR (CDCl3, 400 MHz): δ 1.15 (t, J = 7.0 Hz, 3 H), 1.18 (s, 9 H), 1.21 (s, 9 H), 1.29 (t, J = 7.0 Hz, 3 H), 1.67 (d, J = 7.1 Hz, 3 H), 2.33 (s, 3 H), 3.49 (d, J = 27.6 Hz, 1 H), 3.48 -3.60 (m, 2 H), 3.86 -4.03 (m, 2 H), 4.05 (s, 6 H), 5.44 (q, J = 7.0 Hz, 1 H), 7.13 (d, J = 7.9 Hz, 2 H), 7.70 (s, 2 H), 7.78 (d, J = 8.1 Hz, 2 H); 13 C{ 1 H}-NMR (CDCl3, 75 MHz): δ 15.7 (d, J = 6.6 Hz), 16.0 (d, J = 5.0 Hz),

  , J = 5.5 Hz), 34.8 (d, J = 3.9 Hz), 36.0 (s), 59.8 (d, J = 7.2 Hz), 60.5 (d, J = 7.2 Hz),

  ; 31 P{ 1 H}-NMR (CDCl3, 162 MHz): δ 23.59; HRMS (ESI-TOF) m/z: [M-TsO -] + Calcd for C20H41N3O4P 418.2829; Found 418.2829. imidazol-1-yl)ethoxy)(tert-butyl)amino)-2,2-

4. 2 . 4 . 3 . 4 . 4 . 4 . 5 . 4 . 6 .

 243444546 pKa measurement. pKa values of 2a-b, 3a and 4a-b were measured by monitoring the dependency of 1 H NMR chemical shift at various pH. A 0.01 M solution of alkoxyamine in D2O/CD3OD (1/1) was used. The pH* values, where the asterisk (*) means that the pH was measured in D2O/CD3OD using a pH-meter calibrated with non-deuterated standard solutions, were adjusted with DCl and NaOD to various pH values and 1 H-NMR spectra were recorded on a 400 MHz NMR spectrometer. The chemical shift of proton signals on the imidazole ring was plotted against pH* (see Figure S1-S5 in supplementary data) and a sigmoidal fitting was performed. The obtained pKa* from the curve fitting -which value was determined from pH* -was then corrected to pKa using the following equation pH = 0.929 pH* + 0.42 [11]. Kinetics measurement. The values of homolysis rate constants kd of 2a-c and 4a-c were determined by monitoring the concentration of nitroxide by EPR. A sample tube containing a solution of 1 x 10 -4 M of alkoxyamine in H2O/MeOH (1:1) buffered with Na2CO3/NaHCO3 for pH = 10.3, 50 mM HEPES buffer containing 0.15 mM NaCl for pH 7.4 or HCl/KCl for pH = 2.2 was set in EPR cavity and EPR spectrum were recorded sequentially. The temperature in the cavity was controlled by a temperature controlling unit. The area of the double-integrated signal (I) was plotted against time, and the curve fitting was performed by the exponential fitting to obtain the plateau value (I∞) of the double-integrated signal. Then, the value of ln{(I∞ -I) / I∞} was plotted against time, and this plot was fitted by linear regression to give the kd value. The kd values of 3a-b were determined by monitoring the decay of alkoxyamine by 31 P{ 1 H}-NMR in the presence of TEMPOL (4-hydroxy-2,2,6,6tetramethylpiperidin-N-oxyl) as alkyl radical scavenger to avoid the reverse reaction for the production of the alkoxyamine. A stock solution of 0.02 M of alkoxyamine in H2O/MeOH with 2 equiv. of TEMPOL was prepared. Same buffer solutions as EPR measurement were used for specific pH conditions. The stock solution was sampled into 20 NMR tubes (0.5 mL in each tube). They were sunk in a pre-heated oil bath at reaction temperature, then withdrawn at various time intervals and quenched on a ice-water bath. Then 0.1 mL of DMSO-d6 containing (EtO)3PO (0.002 M) as internal standard (δ = 0 ppm) was added to each tube. P{ 1 H}-NMR signal was recorded on a 400 MHz NMR spectrometer. The value of ln(Int / Int0), where Int is the integral value and Int0 is the integral value at t = 0, was plotted against time and the linear fitting was performed to give the kd value. In general, more than 90% conversion was reached for both EPR and NMR measurements. The kd values were then converted to activation energy (Ea) using Arrhenius equation with the averaged frequency factor A = 2.4.10 14 s -1 . The half-life time t1/2 was calculated from the following equation, t1/2 = ln2 / kd. Cell Cultures. Human glioblastoma cell lines U87-MG an U251-MG were cultured in Dulbecco's modified Eagle's medium (DMEM) and Eagle's minimum essential medium (EMEM), respectively, supplemented with 10% fetal bovine serum (FBS), 100 U/mL of Penicillin and 100 µg/mL of Streptomycin (Pen Strep). EMEM was additionally supplemented with 2 µM of L-glutamine. Cell lines were regularly screened and are free from mycoplasma contamination. Viability assays. Stock solutions of each alkoxyamine were prepared in DMSO at 0.15 M. Cells were treated with a range of concentrations of alkoxyamines (500, 250, 100, 50, 25, 10 and 5 µM) for 72 h. The highest concentration of DMSO of which cells were exposed to was 0.002%. Cells were seeded in 96-well plates (5000 cells/plate) and treated with alkoxyamines. After a 72 h-treatment, cell survival was measured by incubating the glioblastoma cells for 1 h in the medium replenished with 0.5 mg/mL of MTT (3-(4,5dimethylthiazol-2-yl)-2,5-diphenyl tetrazolium bromide). Medium was replaced by DMSO and optical density (OD) was measured by a microplate spectrophotometer at 550 nm. For each dose of each alkoxyamine tested, results were expressed as percentage of survival cells, according to the following equation: ODtreated cells x 100/ODcontrol cells. The IC50 values, i.e.concentrations that inhibit 50% of cell survival, were graphically determined. All experiments were performed at least in triplicate. The viability test in the presence of alkyl scavenger was carried out with an IC50 concentrations of alkoxyamines and an incubation time of 72h. U87-MG and U251-MG cell viability was measured versus equal concentrations of the nontoxic radical scavenger trihydroxyethylrutin (Troxerutin). Images of the treated cells were captured with a camera on an inverted light microscope. LogD7.4 measurements. LogD7.4 values were determined by standard shake-flask procedure using pre-saturated solutions (n-octanol and PBS at pH 7.4). Aqueous phases were then analyzed by RP-HPLC (column: 4.6 × 50 mm, 2.7 µM; solvent: CH3CN (+ 0.1% TFA) / H2O (+ 0.1% TFA), 10:90 to 100/0 in 10 min).

Table 1

 1

	Alkoxyamine	pKa	conditions a	Ea (kJ/mol) b	t1/2 (37 ºC) c	Alkoxyamine	pKa	conditions a	Ea (kJ/mol) b	t1/2 (37 ºC) c
	(RR/SS)-2a	6.09	pH = 10.3	121.6	10 d	(RS/SR)-2a	6.04	pH = 10.3	122.3	13 d
			pH = 2.2	111.1	4 h			pH = 2.2	111.7	5 h
	(RR/SS)-2b	6.28	pH = 10.3	120.0	130 h	(RS/SR)-2b	6.50	pH = 10.3	116.9	39 h
			pH = 2.2	109.7	144 min			pH = 2.2	109.2	119 min
	(RR/SS)-3a	5.62	DMSO	148.9	> 1000 y	(RS/SR)-3a	5.81	DMSO	144.9	> 200 y
			DMSO + TsOH d	142.8	> 100 y			DMSO + TsOH d	140.9	> 50 y
	(RR/SS)-4a	6.11	pH = 10.3	123.0	17 d	(RS/SR)-4a	6.13	pH = 10.3	123.3	19 d
			pH = 2.2	119.4	103 h			pH = 2.2	121.2	8 d
	(RR/SS)-4b	6.76	pH = 10.3	126.4	65 d	(RS/SR)-4b	6.56	pH = 10.3	126.6	70 d
			pH = 2.2	118.2	65 h			pH = 2.2	120.5	7 d

pKa, activation energies (Ea), and half-life time (t1/2) at 37 °C for alkoxyamines (RR/SS)-and (RS/SR)-2a-b, 3a, 4a-b at different pH (see supplementary material for details).

a Solutions were buffered with Na2CO3/NaHCO3 (pH = 10.3) or HCl/KCl (pH = 2.2). b Calculated with experimental kd values using Arrhenius equation and A = 2.4.10 14 s -1 . c Calculated at 37 ºC using re-estimated kd values. d TsOH = para-toluenesulfonic acid.

At basic pH, (RR/SS)-2a and (RS/SR)-2a, as well as (RR/SS)-4a and (RS/SR)-4a, exhibited Ea

Table 2

 2 Activation energies (Ea) and half-life time (t1/2) at 37 °C for alkoxyamines (RR/SS)-

	(RR/SS)-2a	pH = 7.4	118.0	60 h	(RS/SR)-2a	pH = 7.4	118.8	81 h
	(RR/SS)-2b	pH = 7.4	116.3	31 h	(RS/SR)-2b	pH = 7.4	113.5	10 h
	(RR/SS)-2c	pH = 7.4	113.5	10 h	(RS/SR)-2c	pH = 7.4	110.6	204 min
	(RR/SS)-3a	pH = 7.4	-d	-d	(RS/SR)-3a	pH = 7.4	-d	-d
	(RR/SS)-3b	pH = 7.4	135.2	5 y	(RS/SR)-3b	pH = 7.4	135.2	5 y
	(RR/SS)-4a	pH = 7.4	122.1	12 d	(RS/SR)-4a	pH = 7.4	122.3	13 d
	(RR/SS)-4b	pH = 7.4	124.7	33 d	(RS/SR)-4b	pH = 7.4	123.6	22 d
	(RR/SS)-4c	pH = 7.4	115.8	25 h	(RS/SR)-4c	pH = 7.4	116.1	29 h

and (RS/SR)-2a-c, 3a-b, 4a-c at pH = 7.4 (see supplementary material for details).

Alkoxyamine conditions a Ea (kJ/mol) b t1/2 (37 ºC)

c

Alkoxyamine conditions a Ea (kJ/mol) b t1/2 (37 ºC) c a Solutions were buffered with HEPES/NaOH buffer (50 mM, pH = 7.4) with 0.15 mM of NaCl.

Table 3

 3 Molecular

	Alkoxyamine	Molecular weight	logD7.4	Alkoxyamine	Molecular weight	logD7.4
		(g.mol -1 )			(g.mol -1 )
	(RR/SS)-2a	389.5	3.6	(RS/SR)-2a	389.5	4.1
	(RR/SS)-2b	403.5	3.4	(RS/SR)-2b	403.5	3.7
	(RR/SS)-2c	418.5	-0.2	(RS/SR)-2c	418.5	0.5
	(RR/SS)-3a	389.5	3.4	(RS/SR)-3a	389.5	3.4
	(RR/SS)-3b	404.5	-0.1	(RS/SR)-3b	404.5	0.3
	(RR/SS)-4a	389.5	3.7	(RS/SR)-4a	389.5	4.0
	(RR/SS)-4b	403.5	3.5	(RS/SR)-4b	403.5	3.8
	(RR/SS)-4c	418.5	0.8	(RS/SR)-4c	418.5	0.4

weight and measured n-octanol/PBS partition coefficient (logD7.4) of alkoxyamines 2a-c, 3a-b and 4a-c.

  J = 5.5 Hz), 32.7 (s), 35.1 (d, J = 5.5 Hz), 58.6 (d, J = 7.7 Hz), 61.21 (d, J = 6.1 Hz),

	61.22 (s), 70.2 (d, J = 138.1 Hz), 75.7 (s), 120.2 (s), 127.0 (s), 149.2 (s); 31 P{ 1 H}-NMR
	(CDCl3, 162 MHz): δ 25.45; HRMS (ESI-TOF) m/z: [M+H] + Calcd for C19H39N3O4P
	404.2673; Found 404.2672. (RS/SR)-2b: pale yellow oil; Rf = 0.43 (CH2Cl2: acetone = 1: 4);
	1 H-NMR (CDCl3, 400 MHz): δ 1.00 (t, J = 7.1 Hz, 3 H), 1.21 (s, 18 H), 1.26 (t, J = 7.2 Hz, 3
	H), 1.67 (d, J = 6.6 Hz, 3 H), 3.24 -3.41 (m, 2 H), 3.39 (d, J = 25.7 Hz, 1 H), 3.84 (s, 3 H),
	3.85 -3.93 (m, 1 H), 3.93 -4.05 (m, 1 H), 5.36 (q, J = 6.6 Hz, 1 H), 6.73 (d, J = 0.9 Hz, 1 H),
	7.02 (d, J = 0.7 Hz, 1 H); 13 C{ 1 H}-NMR (CDCl3, 75 MHz): δ 15.7 (d, J = 7.2 Hz), 16.0 (d, J
	= 6.1 Hz), 18.5 (s), 27.8 (s), 29.8 (d, J = 6.1 Hz), 32.5 (s), 34.8 (d, J = 5.5 Hz), 58.2 (d, J =
	7.7 Hz), 60.5 (d, J = 6.1 Hz), 61.0 (s), 68.0 (s), 69.9 (d, J = 138.6 Hz), 119.8 (s), 127.5 (s),
	147.3 (s); 31 P{
	4: 1 mixture of diastereoisomers ( 31 P-NMR
	ratio). The crude was separated by flash column chromatography (CH2Cl2: acetone = 85: 15
	to 50: 50) to afford (RR/SS)-2b (1.50 g, 19%) and (RS/SR)-2b (3.11 g, 39%). (RR/SS)-2b: pale
	yellow oil; Rf = 0.33 (CH2Cl2: acetone = 1: 4); 1 H-NMR (CDCl3, 400 MHz): δ 0.84 (s, 9 H),
	1.24 (s, 9 H), 1.31 (t, J = 7.1 Hz, 6 H), 1.69 (d, J = 7.1 Hz, 3 H), 3.32 (d, J = 25.4 Hz, 1 H),
	3.78 (s, 3 H), 3.91 -4.03 (m, 1 H), 4.03 -4.11 (m, 1 H), 4.11 -4.22 (m, 1 H), 4.24 -4.37 (m, 1

H), 5.34 (q, J = 7.0 Hz, 1 H), 6.74 (d, J = 1.0 Hz, 1 H), 6.97 (d, J = 1.0 Hz, 1 H);

13

C{ 1 H}-NMR (CDCl3, 75 MHz): δ 15.8 (d, J = 6.6 Hz), 16.2 (d, J = 6.1 Hz), 19.9 (s), 27.4 (s), 29.5 (d, 1 H}-NMR (CDCl3, 162 MHz): δ 25.06; HRMS (ESI-TOF) m/z: [M+H] + Calcd for C19H39N3O4P 404.2673; Found 404.2671.

  47 g, 19%) and (RS/SR)-3a (1.40 g, 18%). (RR/SS)-3a: pale yellow oil. Rf = 0.23

	(q, J = 6.1 Hz, 1 H), 7.12 (d, J = 7.9 Hz, 2 H), 7.36 (d, J = 1.6 Hz, 1 H), 7.54 (s, 1 H), 7.76
	(d, J = 8.2 Hz, 2 H), 9.96 (s, 1 H); 13 C{ 1 H}-NMR (CDCl3, 75 MHz): δ 15.9 (d, J = 6.6 Hz),
	(CH2Cl2: MeOH = 19: 1); 1 H-NMR (CDCl3, 400 MHz): δ 0.78 (s, 9 H), 1.26 (s, 9 H), 1.33 (t, 16.2 (d, J = 5.5 Hz), 19.9 (s), 20.9 (s), 26.8 (s), 30.3 (d, J = 5.5 Hz), 35.0 (d, J = 4.4 Hz), 36.2
	J = 7.1 Hz, 3 H), 1.34 (t, J = 7.1 Hz, 3 H), 1.78 (d, J = 6.1 Hz, 3 H), 3.39 (d, J = 26.7 Hz, 1 (s), 59.5 (d, J = 7.7 Hz), 61.5 (d, J = 6.6 Hz), 61.6 (s), 69.0 (d, J = 140.3 Hz), 88.9 (s), 120.6
	H), 3.94 -4.05 (m, 1 H), 4.05 -4.17 (m, 2 H), 4.29 (dquint, J = 10.3, 7.3, 7.3, 7.3, 7.3 Hz, 1 (s), 123.7 (s), 125.6 (s), 128.2 (s), 136.9 (s), 138.8 (s), 143.8 (s); 31 P{ 1 H}-NMR (CDCl3, 162
	H), 5.86 (q, J = 6.1 Hz, 1 H), 7.03 (s, 1 H), 7.08 (s, 1 H), 7.72 (s, 1 H); 13 C{ 1 H}-NMR MHz): δ 23.32; HRMS (ESI-TOF) m/z: [M-TsO -] + Calcd for C19H39N3O4P 404.2673; Found
	(CDCl3, 75 MHz): δ 15.9 (d, J = 6.6 Hz), 16.2 (d, J = 6.1 Hz), 19.8 (s), 26.6 (s), 30.5 (d, J = 404.2675.
	6.1 Hz), 35.0 (d, J = 5.0 Hz), 59.1 (d, J = 7.2 Hz), 61.17 (d, J = 7.2 Hz), 61.23 (s), 69.1 (d, J 4.1.8. 1-(1-((tert-Butyl(1-(diethoxyphosphoryl)-2,2-dimethylpropyl)amino)oxy)ethyl)-3-
	= 139.8 Hz), 84.8 (s), 117.0 (s), 128.5 (s), 136.7 (s); 31 P{ 1 H}-NMR (CDCl3, 162 MHz): δ methyl-1H-imidazol-3-ium 4-methylbenzenesulfonate ((RS/SR)-3b). (RS/SR)-3b was
	24.49; HRMS (ESI-TOF) m/z: [M+H] + Calcd for C18H37N3O4P 390.2516; Found 390.2520. synthesized using same procedure as (RR/SS)-2c with (RS/SR)-3a (0.35 g, 0.90 mmol) and
	methyl p-toluenesulfonate (0.94 mL, 1.0 M in THF) to give (RS/SR)-3b (0.10 g, 19%). White
	solid; m.p.: 114-116 o C; 1 H-NMR (CDCl3, 400 MHz): δ 0.85 (br. s, 9 H), 1.21 (br. s, 9 H),
	1.32 (t, J = 7.0 Hz, 6 H), 1.77 (d, J = 6.0 Hz, 3 H), 2.34 (s, 3 H), 3.26 (d, J = 27.1 Hz, 1 H),
	3.89 -4.11 (m, 2 H), 4.03 (br. s, 3 H), 4.12 -4.29 (m, 2 H), 6.29 (br. s, 1 H), 7.14 (d, J = 7.9
	Hz, 2 H), 7.19 (br. s, 1 H), 7.81 (d, J = 8.2 Hz, 2 H), 7.98 (s, 1 H), 10.13 (s, 1 H); 13 C{ 1 H}-
	NMR (CDCl3, 75 MHz): δ 15.8 (d, J = 7.2 Hz), 16.0 (d, J = 6.1 Hz), 20.8 (s), 21.0 (br. s),
	27.4 (s), 29.2 (br. s), 34.8 (br. s), 36.0 (s), 59.4 (br. s), 61.3 (d, J = 7.2 Hz), 63.0 (br. s), 67.6
	(d, J = 135.9 Hz), 89.2 (br. s), 121.5 (br. s), 122.6 (s), 125.5 (s), 128.0 (s), 136.9 (s), 138.4 (s),
	4.1.7. 1-(1-((tert-Butyl(1-(diethoxyphosphoryl)-2,2-dimethylpropyl)amino)oxy)ethyl)-3-143.9 (s); 31 P{ 1 H}-NMR (CDCl3, 162 MHz): δ 24.03; HRMS (ESI-TOF) m/z: [M-TsO -] +
	methyl-1H-imidazol-3-ium 4-methylbenzenesulfonate ((RR/SS)-3b). (RR/SS)-3b was Calcd for C19H39N3O4P 404.2673; Found 404.2676.
	synthesized using same procedure as (RR/SS)-2c with (RR/SS)-3a (0.30 g, 0.77 mmol) and
	methyl p-toluenesulfonate (0.7 mL, 1.0 M in THF) to give (RR/SS)-3b (0.35 g, 79%). Pale
	yellow oil; 1 H-NMR (CDCl3, 400 MHz): δ 0.79 (s, 9 H), 1.19 (s, 9 H), 1.29 (t, J = 7.1 Hz, 3
	H), 1.32 (t, J = 7.1 Hz, 3 H), 1.88 (d, J = 6.1 Hz, 3 H), 2.32 (s, 3 H), 3.35 (d, J = 27.0 Hz, 1
	H), 3.93 -4.14 (m, 3 H), 4.09 (s, 3 H), 4.22 (dquint, J = 10.2, 7.3, 7.3, 7.3, 7.3 Hz, 1 H), 5.97

(RS/SR)-3a: white crystal (recrystallized from EtOAc); m.p.: 68-70 o C; Rf = 0.17 (CH2Cl2: MeOH = 19: 1); 1 H-NMR (DMSO-d6, 500 MHz, 369 K): δ 1.09 (s, 9 H), 1.11 (s, 9 H), 1.14 (t, J = 7.3 Hz, 3 H), 1.26 (t, J = 7.0 Hz, 3 H), 1.67 (d, J = 6.4 Hz, 3 H), 3.37 (d, J = 25.0 Hz, 1 H), 3.60 -4.07 (m, 4 H), 6.10 (q, J = 6.1 Hz, 1 H), 6.86 (s, 1 H), 7.32 (s, 1 H), 7.83 (s, 1 H); 13 C{ 1 H}-NMR (CDCl3, 75 MHz): δ 16.0 (d, J = 6.6 Hz), 16.1 (d, J = 6.1 Hz), 19.4 (s), 27.7 (s), 30.1 (br. s), 35.0 (br. s), 58.9 (br. s), 60.6 (br. s), 61.5 (s), 69.3 (d, J = 140.3 Hz), 83.2 (br. s), 117.1 (br. s), 128.6 (br. s), 137.0 (br. s); 31 P{ 1 H}-NMR (CDCl3, 162 MHz): δ 23.91; HRMS (ESI-TOF) m/z: [M+H] + Calcd for C18H37N3O4P 390.2516; Found 390.2519.

Acknowledgements

The authors were financially supported by Aix-Marseille Université, Centre National de la Recherche Scientifique (CNRS), Institut National de la Santé et de la Recherche Médicale (INSERM) and Assistance Publique des Hôpitaux de Marseille (APHM). The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" program, from the Fondation ARC pour la recherche sur le cancer (nº PJA 20141201886), and from DGOS (SIRIC label) INCa-DGOS-INSERM 6038. This work was also supported by the foundations ReSOP "Réseau d'Onco-hématologie Pédiatrique" and AROU "Association pour la Recherche en Oncologie pédiatrique et Urologique", as well as by the computing facilities of the CRCMM 'Centre Régional de Compétences en Modélisation Moléculaire de Marseille'. The authors thank Pr. Didier Siri for his kind assistance. Marie-Pierre Montero is gratefully acknowledged for her help in the culture of cancer cell lines. Appendix A. Supplementary data

33

Supplementary material to this article (characterization data, kinetic and pKa data, calculations details and x-ray crystallographic data) can be found online at .

Notes and references

1 W. J. Geldenhuys, A. S. Mohammad, C. E. Adkins, P. R. Lockman. Molecular determinants of blood-brain barrier permeation. Ther. Deliv., 2015, 6, 961-971.


Mehdi Bouatou

Rishav Harsh

Cyril Chacon

Yann Girard

Vincent Repain

Amandine Bellec

Sylvie Rousset

Alexander Smogunov

Yannick J Dappe

Jérôme Lagoute
email: jerome.lagoute@u-paris.fr

A Rousset

Black Phosphorus for directed molecular assembly with weak electronic coupling

Keywords: scanning tunneling microscopy, black phosphorus, molecules, density functional theory

The combination of two-dimensional materials with organic molecules offers the possibility to obtain low dimensional hybrid materials with tailored properties. Black phosphorus is a monoelemental semiconducting material with a non-planar puckered atomic structure. Here, we study porphyrin molecules on black phosphorus and show that they self-assemble with an epitaxial relationship indicating that the molecule-surface interaction largely dominates the molecule-molecule interaction. The atomic structure of black phosphorus is found to be at the origin of the substrate driven self-assembly. Despite this strong interaction with molecules, the electronic coupling is found to be weak allowing the molecules to maintain the properties of their gas phase, which is the usual behaviour for van der Waals materials. Therefore, the combination of the peculiar puckered structure with the van der Waals nature of black phosphorus provides this material with the ability to interact strongly enough with adsorbed molecules to drive their assembly but weakly enough to keep their electronic properties intact.

Introduction

The combination of organic molecules with two-dimensional (2D) materials provides a rich playground for the realization of hybrid materials with improved performances beyond the pristine materials. The molecule-substrate interface is the place where several physical phenomena occur such as charge transfer [1], gap renormalization [2], self-assembly [3], redox reactions [4]. Mastering these interactions would open the route toward interfacial engineering with various potential applications. Molecule-graphene hybrids have been used to create a bandgap in graphene with potential applications in graphene based electronic devices [5]. The functionalization of graphene with metalloporphyrins has been demonstrated to improve the performances of graphene based photodetectors and chemical sensors [6]. On black phosphorus (BP), non-covalent functionalization has been achieved with expected applications in electronics, optoelectronics, energy storage, sensors, and fillers for composite reinforcement [7]. Non-covalently or covalently bonded molecules on BP have been used for surface passivation [8,9,10,11] together with the formation of ohmic contact and charge transfer that are ideal for BP based transport devices [12]. On 2D van der Waals materials, the molecule-surface interaction is generally weak [13], which has two major consequences. Firstly, the weak electronic coupling allows the molecules to maintain the electronic properties of their gas phase [13]. Secondly, molecular assembly on these materials can occur independently of the structure of the 2D material, except when the properties of the 2D materials are modified by the interaction with a substrate which often leads to the formation of moiré patterns [14,15]. Indeed, molecular assembly on surfaces depends on the interplay between molecule-molecule and moleculesurface interaction [3,16,13]. On non-interacting surfaces, the assembly is driven by intermolecular interactions which determine the molecular lattice parameters [17,18,19]. On interacting surfaces, the molecule-substrate interaction dominates and the substrate governs the assembly [20,21,22,23]. In principle, on weak interacting surfaces, the forces driving the self-assembly and the electronic coupling are not directly linked since the former involves electrostatic and van der Waals interaction and the later involves hopping integral between molecules and surfaces wavefunctions. However, the case of interacting surfaces, where the substrate drives the molecular assembly is generally found on metallic surfaces, where the electronic coupling is non negligible. On van der Waals materials, the molecule-substrate interaction force is weak, and the electronic interaction is weak. An intermediate regime where a weak electronic coupling would combine with a strong adsorption force would allow to govern the molecular assembly by the substrate geometry while preserving their electronic properties. So far, interfacial effects on molecule/2D hybrids have been explored at the atomic scale with scanning tunneling microscopy (STM) on various substrates including graphene, hexagonal boron nitride and transition metal dichalcogenide [24,13]. However, little is known for molecules on black phosphorus (BP). In contrast with other 2D materials, black phosphorus is a van der Waals material with a non-planar puckered structure. Its physical properties have been thoroughly investigated at the atomic scale by STM [25,26] revealing anisotropic structural and electronic properties. Here, we study the interaction of iron tetraphenylporphyrin (FeTPP) molecules with BP. We show that the puckered nature of the atomic structure of BP provides to this material the property to drive the self-assembly process of organic molecules leading to a commensurable 3 × 4 molecular lattice. We evidence the weak electronic interaction between FeTPP and BP by comparing the properties of FeTPP/BP with FeTPP/graphene, that serves as a reference for non-interacting 2D material, and combining STM and density functional theory (DFT) calculations. Our results show that the van der Waals nature of BP combined with its non-planar structure gives to BP the property to have a strong interaction with molecules that drives the molecular assembly, while preserving the gas phase electronic properties of the molecules.

Results and discussion

A monolayer of FeTPP on BP forms a rectangular lattice as shown in Figure 1a. There are two types of molecules characterized by a bright or dark center. The molecules with a bright center correspond to FeTPP with a chlorine bonded to the Fe atom forming a FeTPP-Cl molecule [27,28]. Those with a dark center correspond to FeTPP. We focus our study on FeTPP only, since the spectroscopic measurements on FeTPP-Cl often induce removal of the Cl ion and tip modification [29]. The unit cell vectors of the molecular lattice have lengths 1.25 ± 0.01 nm along the [100] direction and 1.31 ± 0.01 nm along the [010] direction (see more examples in supporting information Figure S1) which are smaller than the intermolecular distance obtained in FeTPP layers (and other TPP molecules) on other surfaces (see supporting information Figure S2). Indeed, on Au(111) FeTPP arrange in a rectangular lattice with a periodicity of 1.40 ± 0.05 nm and 1.55 ± 0.05 nm [27], on Ag(111) it forms a square lattice with a period of 1.43 ± 0.01 nm [30], on graphene [28] we have found a square lattice with a periodicity of 1.4 ± 0.02 nm (see supporting information Figure S3). Generally, tetraphenylporphyrin molecules adopt a square lattice with a period of 1.40 nm irrespective of the presence and nature of a central metal atom [31,18]. Here the distance between adjacent molecules is shorter than on metallic substrates or graphene, which is due to the interaction with BP that forces the molecules to follow the BP atomic lattice. As a consequence, the density of TPP molecules on BP is increased by 20% with respect to TPP on other substrates. This could be advantageous in TPP based data storage systems such as H 2 TPP on surfaces. Indeed it has been shown recently that the control of tautomerisation in self-assembled H 2 TPP molecules on graphene allows to store binary information on each molecule [32]. The increase of molecular density on BP can be used to increase the storage density on such systems. It can be noted that the molecular density found here is also higher than tin phtalocyanine (SnPc) that have been used as molecular lattice for data storage on graphene even though single Pc molecules have a smaller size than TPP derivatives [33]. Moreover, due to the epitaxial relationship, TPP derivatives on BP form an addressable system where the position of the molecules is predetermined by the atomic lattice of the substrate. Another implication of the increased molecular density on BP is the potential gain of efficiency in hybrids materials involving TPP derivatives for catalysis. FeTPP is an efficient catalyst for O 2 reduction [34] and CO 2 reduction [35] reactions. The immobilization of porphyrin on 2D materials has been shown to improve the catalytic activity of porphyrin compounds for CO 2 reduction [36]. The increase of in-plane density may further improve the performances of such hybrid materials. This comes from the interaction with the BP substrate which forces the molecules to follow the BP atomic lattice.

In order to understand the role of the BP lattice on the molecular assembly, one has to determine the relationship between the molecules and the underlying BP atomic lattices. Usually, this information is not accessible by STM since the substrate underneath the molecules cannot be visualized. Here, to be able to visualize the substrate below the molecular layer, we have removed all the molecules using molecular manipulation on large scale with the STM tip. We have scanned the surface several times with a reduced tip-sample distance in order to sweep away the molecules. After several scanning with various setpoint parameters, we were able to recover a good tip state allowing to visualize the atomic lattice of the BP on the same area. This procedure is similar to the one that we have used in our previous works on graphene [37,28,4]. The image of the BP surface shows the phosphorus zigzag rows (Figure 1b). Indeed, a single BP layer has a puckered structure made of two planes of phosphorus atoms. In the following we will refer to top and bottom P atoms in a BP layer depending on whether the atoms are in the top (facing the vacuum side) or bottom (facing the bulk side) plane. The atomically resolved STM image corresponds to the top P atoms that form zigzag lines. The measured BP unit cell vectors length are | a| = 4.17 ± 0.1 Å along the armchair direction [100] and b = 3.27 ± 0.1 Å along the zigzag direction [010]. Consequently, the FeTPP layer remarkably matches the BP atomic lattice with a ratio of 3 along the armchair direction and 4 along the zigzag direction, leading to a 3 × 4 pattern. This epitaxial relation is further evidenced in the FFT image of Figure 1c that reveals the periodicity of BP and FeTPP. The FFT was obtained from a combined image made of the addition of the images of Figures 1a and1b. This procedure allows to visualize in one single image the FFT of the BP and the FeTPP lattices and to highlight the 3 × 4 pattern. This epitaxial relationship indicates that the molecule-surface interaction is the dominating force that stabilizes the molecular lattice. The relaxed structure of FeTPP on BP from our ab initio calculations reveals that the molecules adsorb on the surface with the Fe atom located above the topmost zigzag row of BP as shown in Figure 1d. The molecular axis represented in Figure 1d (connecting two opposite phenyl groups) is rotated by an angle of 35.5 • with respect to the molecular lattice vector which is aligned with the zigzag direction. This azimuthal angle corresponds to an intermolecular interaction of "T-type" consisting of an edge to face arrangement of the phenyl groups [31]. In order to understand the driving force for the molecular self-assembly on BP we have performed ab initio calculations including van der Waals forces. First, we have determined the adsorption energy of one FeTPP molecule on the BP surface. To this end, we have considered a 3×4 unit cell of BP composed of 3 atomic layers. The FeTPP molecule has been placed on top of the BP slab, and the full system has been optimized until the forces went below 0.1 eV/ Å using the localized-orbitals basis set DFT code Fireball [38]. From the optimized configuration, we have then calculated the interaction energy between the FeTPP molecule and the BP (defined by the difference in total energies E F eT P P +BP -E F eT P P -E BP ) at different heights from the BP surface. Van der Waals interactions have been taken into account through the use of a specific formalism, the LCAO-S 2 +vdW approach [39], previously used for similar systems [40,41]. The calculated interaction energies are presented in Figure 2a. We compare this interaction energy with the case of FeTPP on graphene that we take as a reference for weakly coupled FeTPP [28]. The optimized equilibrum configuration on BP is presented in Figure 2b. As it can be clearly seen, the adsorption energy on BP is much higher (about 2.4 eV) than the one obtained on graphene (about 1.5 eV). Note also that the most stable position of the molecule on BP is obtained for the molecule center on a zigzag row of top P atoms. Secondly, we have determined the interaction energy of a molecule in the network of its neighboring molecules. Using the same approach as for the interaction with the BP surface, we have considered a set of nine molecules isolated in a plane, and we have calculated the interaction of the central molecule with its neighbors as a function of the intermolecular distance. The result is represented in Figure 2c, and the corresponding molecular network at equilibrium is represented in Figure 2d. The obtained interaction energy is about 0.78 eV per molecule, namely three times less than the interaction of the molecule with the BP surface. Also the intermolecular distance in a free standing layer obtained in this calculations is 13.85 Å in one direction and 14.35 Å in the other one, in both cases much bigger than the lattice parameters determined experimentally for the molecules on BP. Therefore, the molecular network is constrained upon adsorption on the BP surface. This constraint can be roughly estimated energetically at around 2.5 eV loss in the interaction energy between molecules, when extrapo-lating the curve in Figure 2c. This could appear to be huge, but it has to be compared to the interaction energy with the substrate presented in Figure 2a which is about 2 eV. Moreover, we have to stress that the results presented in Figure 2c arise from a static calculation, starting from fully optimized isolated molecules, arranged in a network and displaced rigidly without further optimization. Obviously, the very high repulsion found in Figure 2c would be attenuated through a further optimization, leading to a more important substrate-molecule interaction than the molecule-molecule repulsion. In other words, these findings fully confirm the fact that the BP substrate stabilizes the molecular network in a constrained phase. The surface-driven molecular organization is due to the unique puckered atomic structure of BP. Indeed, along the armchair direction, the van der Waals interaction results in high energy gain when the molecule locates above the top P atoms because the molecule-phosphorus distance is much shorter as compared to the bottom P atoms. This allows to understand that along the armchair direction, the molecules match the BP atomic lattice. Along the zigzag direction, a first intuitive argument would be that all the P atoms are at the same height and therefore the molecules should not be constrained in that direction. In order to gain a deeper insight into the driving force of the self-organization, we plot the charge redistribution between BP and FeTPP in Figure 2e,f. In Fig 2e the isosurface of the charge density difference shows that above BP the charge redistribution is dominated by the top P atoms. This charge redistribution is a consequence of the puckered structure of BP that leads to a reduced interaction between bottom P atoms and FeTPP as compared to the top P atoms. This is in line with the argument given above that the van der Waals interaction is weaker for the bottom P atoms. Strikingly, in contrary to the first intuitive expectation, a strong corrugation in the charge redistribution is also present along the zigzag direction (Figure 2f). This indicates that the molecule 'see' a non uniform potential also along the zigzag line. This is a consequence of the atomic structure of BP, not only the buckling that occurs in the armchair direction, but also the large interatomic distance between the top P atoms. Indeed, in the [010] direction, the P atoms are separated by a distance of 3.3 Å which is larger than the standard P-P bond-length. As a consequence a strong variation of the potential occurs in that direction and the van der Waals interaction is non uniform leading to a constrain in the molecular assembly in the zigzag direction. Therefore, the constrained self-assembly that occurs on BP is a consequence of the geometry of the atomic lattice that leads to strong periodic variations of the van der Waals interaction between the molecule and the surface along the armchair and zigzag direction. On other 2D materials such as graphene, the shorter interatomic distance prevents such strong variations of the force interaction and charge redistribution leading to a smooth potential variations on the surface and to molecular assembly driven by intermolecular forces which is opposite to the molecule-surface dominating interaction on BP. The BP appears to be an interacting substrate that drives the organization of adsorbed molecules. Therefore a strong electronic interaction between BP and FeTPP is expected. Surprisingly this is not the case. The first information on the electronic properties of FeTPP on BP is obtained by visualizing the HOMO and LUMO states from STM topographical images (Figure 3). The central part of the molecules appears as a protrusion at negative bias (-1 V, Figure 3a,c) and as a depletion at positive bias (+2 V, Figure 3b,d). Note that at negative bias both the FeTPP-Cl molecules and FeTPP appear with a protrusion at the center and that FeTPP appear brighter than FeTPP-Cl unlike what is observed at positive bias. The images of FeTPP on BP are very similar to the ones of FeTPP on graphene [28] at -1 V and +2 V. On graphene the image at -1 V arises from the HOMO state localized on the Fe atom, and the image at +2 V arises from the LUMO and higher unoccupied states localized on the organic part. The same trend is observed here, on BP. A close look at the images in Figure 3c,d reveals four lobes at the periphery of the molecule corresponding to the phenyl groups (as indicated by the molecular model superimposed to the experimental images). The corresponding azimuthal angle θ with the [010] axis of BP measured in Figure 3c is 31±4 • , close to the value obtained in our calculation. From these topographical images we conclude that the HOMO state is localized on the iron atom, and that the azimuthal angle is in agreement with the calculations, which corresponds to a T-type interaction between molecules. More insight is given by the spectroscopy measured at different positions above a FeTPP molecule. In Figure 4, we show the spectra measured above the Fe atom, a pyrrole group, and a phenyl group. Two peaks are visible at negative bias corresponding to the HOMO-1 at -1.38 V, and the HOMO at -1 V.

At positive bias, a broad peak appears for which 2 contributions can be distinguished at 1.95 V and 2.2 V. The measured HOMO-LUMO gap is 2.95 V. The peak of the HOMO state dominates the spectrum at the center of the molecule, revealing that the HOMO state is localized on the Fe atom. The HOMO-1 and the peaks at positive bias are localized on the pyrrole and phenyl groups. The localization of the states can be further evidenced from dI/dV maps which are shown in Figure 4d. These electronic properties are very similar to those of FeTPP on graphene [28] (see supporting information Figure S4) where a HOMO-LUMO gap of 2.9 V was measured and the localization of the molecule electronic states was similar. A noticeable difference is found at positive bias where 3 sharp peaks could be clearly distinguished on graphene around the LUMO state, whereas a broad peak is observed here and only two contributions can be distinguished. This discrepancy comes from a difference in the electronic properties of graphene and BP. While the local density of states (LDOS) of graphene increases linearly with the energy and symmetrically around the Dirac point, the LDOS of BP exhibits a bandgap but also a strong asymmetry with a strong increase of the measured LDOS at positive sample bias voltage as shown in Figure 4c. Around the Fermi level, the spectrum shows a gap of 0.3 V with a valence band edge close to the Fermi level indicating that the BP is p-doped. The measured gap corresponds to the expected value for bulk BP and the spectroscopy is in agreement with previous STM reports [25,26]. At larger bias, we observe the presence of a peak at 1.5 V. The asymmetry of the dI/dV spectrum with a large value at positive bias explains that the molecular states around the LUMO undergo a large broadening as compared to the states at negative bias. This explains why on BP the molecular spectrum around the LUMO state is not as resolved as on graphene.

To interpret the experimental spectra we have carried out DFT calculations using plane wave electronic structure package Quantum ESPRESSO (QE) within local density approximation (in Perdew-Zunger parametrization) for exchange-correlation potential. We have applied a Hubbard-like correction of U = 5 eV to Fe d orbitals, as in our previous paper for FeTPP on graphene and Au(111) [28], in order to compare to those results. The calculations for FeTPP on three-layers BP slab reveal that the lowest energy state has a spin state S=1 with a (d xy ) 2 (d xz ) 1 (d yz ) 1 (d z 2 ) 2 configuration as demonstrated in Figure 5 where total, molecular and Fe d-orbital resolved DOS are presented. This state is generally found to be the ground state of isolated FeTPP [42,43,44,45]. On BP, this state is lower in energy by about 80 meV compared to other S=1 states with configurations (d xy ) 2 (d xz ) 2 (d yz ) 1 (d z 2 ) 1 or (d xy ) 2 (d xz ) 1 (d yz ) 2 (d z 2 ) 1 , which is similar to the case of graphene [28] (but different from Au(111)) indicating relatively weak electronic hybridization of molecular states with BP substrate. One can observe, in particular, that HOMO, located at about 1 eV below the Fermi energy, is heavily dominated by the mentioned above Fe (spin down) d z 2 orbital (shown in green) which extends largely in the vacuum and can be therefore readily probed by STM, in agreement with experimental observations. Unocuppied molecular states are, on the contrary, localized on molecular skeleton (which also agrees with the experimental STS spectra) and are non spin-polarized. The present results are very close to those for FeTPP on graphene and to the fundamental state of a free FeTPP molecule [28]. These electronic properties, together with the self-assembly described above show that BP is a material that interacts strongly enough with molecules to govern their assembly, but weakly enough to preserve their intrinsic electronic properties. In that respect, BP is a unique material that can combine these 2 properties which usually exclude each other.

Conclusion

In conclusion, we have studied the interaction of black phosphorus with FeTPP molecules. The molecules self-assemble on the surface following the periodicity of BP. This assembly is due to the molecule-surface interaction that dominates the molecule-molecule interaction. The unique structure of BP combining zigzag and armchair rows leads to an energy gain of the system due to van der Waals interaction that dominate the energy loss induced by the reduced intermolecular distances. As a consequence, the molecules adopt a dense lattice to fit with the periodicity of BP. The elecronic spectra of the adsorbed molecules reveal a S=1 groundstate similar to the one previously observed on graphene and to the free molecule. These results show that BP allows to drive molecular self-assembly while preserving the electronic prop-erties of the molecules. This opens new avenues for functional surfaces based on 2D materials.

Experimental Section

Black phosphorus samples were purchased from HQ Graphene. The samples were cleaved under ultrahigh vacuum conditions. All STM measurements were performed with a low-temperature STM apparatus (Omicron) working at 4.6 K at a pressure lower than 1 × 10 -10 mbar. The dI/dV spectra were acquired using a lock-in detector at a frequency of ca. 823 Hz and a modulation amplitude of 35 mV. The measurements were performed with an electrochemically etched tungsten tip. Before measuring on graphene, the tip was calibrated on a Au(111) surface until it showed the Shockley surface state feature in the spectroscopic measurements.

Table of Contents

Iron tetraphenylporphyrin molecules on black phosphorus are studied by scanning tunneling microscopy. The epitaxial relation between molecules and black phorphorus, and the local spectroscopy reveals that the interaction of black phosphorus with molecules is strong enough to drive their assembly while the electronic coupling is weak enough to allow the molecules to maintain the electronic properties of their gas phase.

Figure 1 :

 1 Figure 1: Assembly of FeTPP on BP. (a) STM image (2 V, 10 pA) of self-assembled FeTPP molecules adsorbed on BP. The black arrows indicate the crystallographic directions of BP, the yellow arrows and square indicate the unit cell of the FeTPP lattice. (b) STM image (-0.4 V, 200 pA) of the same area as in (a) after removing the molecules with the STM tip. The inset is a zoom of atomically resolved image with the BP unit cell indicated. (c) FFT of the superposition of the STM image of the FeTPP molecules with the image of the BP atomic structure. The peaks corresponding to the FeTPP lattice are marked with red circle, those of the BP atomic lattice are marked with black circles. Intermediate spots along the two main axis are indicated with white markers. (d) Relaxed structure of FeTPP on a BP monolayer. Color code: Black for phosphorus, yellow for carbon, blue for nitrogen, red for Fe, white for hydrogen.

Figure 2 :

 2 Figure 2: (a) Interaction energy of FeTPP on black phosphorus as a function of the molecule-surface distance, compared to the interaction energy for FeTPP on graphene. (b) Side view of the optimized structure of FeTPP on BP. (c) Interaction energy of one FeTPP molecule in the middle of the self-assembled network of its neighboring molecules, as a function of the Fe-Fe distance between neighboring molecules showing a minimum for a Fe-Fe distance of 14.35 Å. (d) Structure of the FeTPP layer considered for the interaction energy calculated in (c). (e,f) Side view of the charge density difference (isovalue 0.01 e).

Figure 3 :

 3 Figure 3: Bias voltage dependent STM images of FeTPP. STM images of FeTPP on BP measured at (a) -1 V, 10 pA and (b) 2 V, 10 pA. (c,d) Zoomed images of (a,b) on FeTPP without Cl atom with a superimposed model of FeTPP as a guide to the eye.

Figure 4 :

 4 Figure 4: Electronic properties of FeTPP on BP. (a) STM image (2 V, 40 pA) of a FeTPP molecule with a model that serves as a guide to the eye. (b) dI/dV spectra measured at the locations indicated by crosses in (a). (c) dI/dV spectrum of BP. The inset is a close view in a small voltage range around the Fermi level showing the bandgap of BP. (d) Topography (2 V, 60 pA) and dI/dV maps of four FeTPP molecules on BP (image size 2.8 nm × 2.8 nm).

Figure 5 :

 5 Figure 5: Spin-resolved DOS for FeTPP on three-layer slab of BP. Top: total DOS (black), as well as DOS projected on FeTPP (red) or on Fe (blue). Bottom: PDOS (projected DOS) on different Fe d-orbitals. The right panels show the relaxed geometry of the FeTPP on BP.

Acknowledgements

This project received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement No. 766726. DFT calculations have been performed using HPC resources from GENCI-CINES (Grant A0100910407).

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.