Siffeen Zehra

Huzaifa Yasir Khan

Thierry Roisnel

Sartaj Tabassum

Farukh Arjmand
email: farukh.arjmand18@gmail.com

Structural insights into interactions of new polymeric (µ-oxo) bridged Cu(II) complexes of taurine with yeast tRNA by spectroscopic and computational approaches and its application towards chemoresistant cancer lines

Keywords: Binding modes, Cu(II) complexes, Drug-RNA interaction, EtBr displacement, cytotoxicity

diffusion de documents scientifiques de 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

Medicinal chemists globally, are in search of ‗robust' efficacious safe chemotherapeutic drugs for cancers that can overcome the existing issues of currently used drugs viz., acquired resistance, off-target toxicity and severe side effects [1][2][START_REF] Munteanu | Metal Anticancer Compounds[END_REF][4]. Many drug design strategies have emerged over a period of time, as most of the existing anticancer drug regimens have failed at one or other stages of R& D or clinical trials that puts huge economic toll on countries and also wasted considerable span of time [5][6][7][8]. Amongst the innovative strategies is the use of RNAbased therapeutics that have recently gained much prominence in drug design for the treatment of chronic diseases such as developmental health issues {metabolic, immunological and neurological}and cancers [9][10][11]. RNAs fold in myriad secondary structural conformations (Hairpin, bulges, loops pseudo knots, etc.) that provide unique binding pockets for selective and specific recognition of small drug molecules [START_REF] Derose | 8.39 -Metal Ion Interactions With DNA, RNA, and Nucleic Acid Enzymes[END_REF][START_REF] Arjmand | [END_REF]. In this context, several intercalative cations, antitumor aminoglycoside antibiotics, amino acids and peptides, inorganic compounds that include Schiff base complex drug formulations derived from many bioactive ligands have been explored and their interaction with the therapeutic targets, DNA and sister cousin, RNAs was vigorously undertaken [14]. Transition metal-based drug entities are unique in geometric shapes [15] and exhibit specific preferential disposition at the biological spaces in vivo, interact readily by both covalent as well as non-covalent interactions [16]. As a result, inorganic drug candidates may be able to modify or tune the conformation of RNA binding sites depending on the metal involved or the bioactive ligand scaffolds, and they may do so with exceptionally specific molecular selectivity and recognition [17].

ACCEPTED MANUSCRIPT / CLEAN COPY

Literature reports reveal that structure-based design motifs are pharmaco-kinetically more druglikely in terms of drug stability, delivery and specificity towards therapeutic targets at the molecular levels [18]. Nevertheless, one of the important identified hurdles for the failure of drugs is poor solubility in water, which limits the evaluation of their pharmacokinetic parameters from in vitro to pre-clinical in vivo research [19]. Literature reports reveal that drugs with low water solubility have a low oral bioavailability that show varied clinical response. To address the solubility issue, medicinal chemists have utilized water soluble analogues of such therapeutics [20] obtained by modifying the core ligand framework by incorporating sulfonic functional groups such as taurine ligands. The sulfur-containing analogue of β-amino acid, taurine (2aminoethane sulfonic acid) is the ubiquitous free amino acid discovered in the mammalian tissues [21]. It is involved in many physiological processes such as cell membrane stabilization, regulating intracellular calcium levels, cell growth and differentiation [22,23]. Taurine and its derivative compounds exhibit remarkable anticancer effects that could be manifested by the following i) suppression of cell proliferation ii) downregulating the angiogenesis process and iii) enhancing the tumor cell apoptosis by scavenging upon the reactive oxygen species [24]. Recent studies have suggested that the intracellular taurine levels could predict the formation of certain tumors and serve as novel biomarkers to diagnose cancers like breast and bladder cancer at an early stage. Therefore, taurine-as a pharmaceutical auxiliary scaffold can be added for better therapeutic response in many classes of anticancer agents [25]. Earlier reports from Schiff bases having taurine moiety exhibit myriad coordination modes while a ternary complex interacts by characteristic aromatic-ring stacking intercalation mode that could stabilize the double-helical structure of therapeutic target, DNA. H. Liang et al. reported a Cu(II) complex, [Cu(saltau(phen)]1.5H 2 O possessing salicylaldehyde-taurine Schiff base and 1,10-Phenanthroline and ACCEPTED MANUSCRIPT / CLEAN COPY studied its DNA interaction behaviour. They demonstrated a dose-dependent antitumor effect on cell cycle and apoptosis in MGC-803 cells [26]. However, there are fewer reports of taurine Schiff base complexes targeted to tRNA that show exceptional anticancer profile. A conjugate is formed by taurine with a uridine residue of tRNA Leu(UUR) in normal mitochondria. The conjugation reaction can improve the interaction of the tRNA's AAU anticodon as the modified uridine residue is situated in the anticodon's wobble location. [27].

Copper complexes have appeared as the most prominent potent class of drugs having anticancer properties [28][29][30] and a series of copper compounds casiopeínas® have entered phase I clinical trials [31], however, they exhibit very low aqueous solubility. Interestingly, copper complexes exhibit a good cytotoxic response against different cancer cell strains and are capable of inhibiting cancerous cells' growth in a dose-dependent manner at a low micro-molar IC 50 values [32]. Copper-based anti-cancer agents are known to kill cancer cells by various mechanisms such as generation of reactive oxygen species (ROS), proteasome inhibition, DNA damage via apoptosis etc [33][34][35][36].

The combinational cocktail targeted therapy incorporating multiple bioactive ligands in a single entity such as electron-rich heterocyclic imidazole as a co-ligand with Schiff base anchoring fragment is a unique approach that has proved instrumental for imparting enhanced cytotoxic activity with reduced toxicities that may necessitate a reduction in dosage of drugs. Previous literature reports revealed that imidazole-based compounds are proven as potent anti-cancer, anti-inflammatory, and antiviral agents [37,38]. Imidazole is an aromatic moiety having both polar and ionizable property that is known to improve pharmacokinetic properties, enhance the solubility and bioavailability of therapeutic drug candidates [39]. The synthesized polymeric copper(II) taurine Schiff base appended with imidazole ligand 1 and 2 drug entities reported, ACCEPTED MANUSCRIPT / CLEAN COPY herein, are exceptionally unique in design, structurally well characterized, specifically targeted to tRNA, fairly soluble in water and capable of inducing ROS mediated cell death in cancer cells.

Experimental

Synthesis of Complexes

To the stirred methanolic solution of 2-Aminoethane sulfonic acid (0.125 g, 1 mmol) was added 2-Hydroxy-3-methoxybenzaldehyde (0.152 g, 1 mmol)/ 3-Hydroxy-2-naphthaldehyde (0.172 g, 1 mmol) and the reaction mixture was kept on reflux for 3 h to yield a yellow-colored solution.

The completion of the reaction was monitored by using TLC and subsequently Cu(CH 3 COO 2 ).H 2 O (0.199g, 1 mmol) was added which resulted in green colored solution. After half an hour a methanolic solution of imidazole (0.068 g, 1 mmol) was added to the reaction. The reaction mixture was filtered and kept aside at room temperature. Suitable green, thin, needle like crystals were obtained from mother liquor after 3-4 days.

Results and discussion

Two polymeric imidazole appended water-soluble Cu(II) complexes 1 and 2 have been prepared as shown in Scheme 1. Light green, thin, needle-like transparent crystal of complex 1 was obtained from the mother liquor, while in the case of complex 2, the compound was crystallized ACCEPTED MANUSCRIPT / CLEAN COPY in a solvent system containing methanol and water (1:1). The spectroscopic techniques (FT-IR, UV-vis and EPR), ESI-MS as well as single X-ray crystallographic studies were utilized to establish the molecular structure of complexes 1 and 2. The complexes were discovered to be airstable and found to be soluble in water, MeOH, EtOH, Dimethylsulfoxide (DMSO) as well as Dimethylformamide (DMF).

ACCEPTED MANUSCRIPT / CLEAN COPY Scheme 1. Synthetic route to prepare the complexes 1 and 2.

Spectroscopic characterization

IR Spectral studies

The complexes 1 and 2 showed the presence of vibration bands associated to ν(-C=N) stretching of the azomethine group centered at 1622 and 1618 cm -1 . The characteristic vibrational frequencies for the stretching modes of ʋ(S=O) group were observed at 1243 and 1139 cm -1 for 1 while at 1236 and 1136 cm -1 for 2 [40]. The absorption bands for the ν (C-N) stretching vibrations corresponding to an imidazole ligand were observed in the range of 1474-1434 cm -1 .

Far-IR spectra of complexes 1 and 2 exhibited peaks in the range of 588-543 cm -1 and 499-486 cm -1 attributed to ν(Cu-N) and ν(Cu-O) stretching vibrations, respectively.

Electronic Spectral studies

The absorption spectra of 1 and 2 were recorded in conc. range 200-700 nm in methanol. The UV-vis spectrum showed bands at 263 and 266 nm in the UV region that were attributed to electronic transitions (π-π*) of the aromatic chromophore while LMCT bands (ligand-to-metal charge transfer) appeared at 347 and 354 nm. These complexes 1 and 2 also displayed lower energy d-d bands due to the 2 Eg → 2 T 2g transition centered at 678 and 680 nm in the visible region. The electronic spectral data indicated that the complexes formed distorted square pyramidal geometry around Cu(II) well-corroborated with single X-ray crystal diffraction studies [41].

EPR Spectral studies

The EPR spectra of complexes 1 and 2 in the polycrystalline form were acquired at room temperature (Fig. S1) to investigate the coordination geometry of the complexes. Complex 1 showed a single broad signal with g∥ at 2.34 and g⊥ at 1.92, whereas complex 2 had g∥ signals at 2.28 and g⊥ at 1.90, which implicated a penta-coordinated geometry around the Cu(II) metal ACCEPTED MANUSCRIPT / CLEAN COPY centers. The penta coordinated complexes with g∥ value greater than g⊥ suggested a distorted square pyramidal environment, ruling out the behaviour of complexes possessing a trigonal bipyramidal geometry [42].

Mass Spectral studies

The ESI-MS spectra of complexes 1 and 2 were recorded, which shows positive ion peaks [M + ] at m/z values corresponding to the molecular formulas of the complexes (Figure S2). The

Single X-ray crystal structure description

The complex 1 crystallizes in the trigonal crystal system with the R-3 space group, according to single X-ray crystal diffraction experiments. The copper metal center adopted a distorted pentacoordinated square pyramidal geometry via coordinating through the N 2 O 3 donor set. The coordination sphere consists of -ONO donor set from the tridentate Schiff base ligand, one nitrogen atom from the imidazole and one oxygen atom from the sulfonate group of the neighboring imine ligand occupying the axial position and acting as a -bridge‖ between the two copper centers arranged perpendicularly (Fig. 1). The oxygen atom at the apex of the pyramid ACCEPTED MANUSCRIPT / CLEAN COPY

Density functional theory computational studies

Density Functional Theory (DFT) computational studies were used to optimize the electronic structure and calculate various geometrical parameters of 1 and 2 using ORCA 3.0.1 program suite. The initial geometries which were subjected to optimization were taken from the X-ray crystallographic data of complexes 1 and 2 (Fig. S3). The associated geometrical parameters, including bond angles as well as bond lengths are determined from the optimized structures and were found to correlate with the crystallographic data as presented in Tables S4 andS5.

However, slight differences in the values resulted from the optimization studies in the gas phase, whereas the experimental data were collected in the solid phase.

The frontier molecular orbitals of a complex were related to optical and absorption properties, within the delocalized orbitals (Fig. 5). The energy for the LUMO and HOMO orbitals were observed in negative values, which demonstrated that the complexes are thermodynamically stable in nature. Further, the orbital analysis for complexes 1 and 2 shows that the E LUMO values were -1.770 and -2.121 eV, whereas E HOMO values were found to be -5.568 and -5.961 eV, respectively. The ΔE values resulted from the DFT calculations were found as -3.798 eV and -3.840 eV for 1 and 2. These values determined that complex 2 had a higher stability in comparison to 1. In conclusion, higher stability of complex 2 and the lower energy band gap values compared to 1 were observed in concordance with the findings observed from various in vitro assays.

Hirshfeld surface analysis

ACCEPTED MANUSCRIPT / CLEAN COPY

The Hirshfeld surface analysis was performed for complexes 1 and 2 as shown in Fig. S4, mapped over d norm , shape index, and curvedness. The d norm surface displayed circular red color depressions which arise from the O----H interactions within both the complexes. The shape index surface demonstrated the red and blue conjoined triangles in the shape of -bow-tie‖ patterns which were attributed to shape the π-π stacking interactions in the crystal structures of the complexes [43].

The two-dimensional fingerprint plots were generated to make it easy to quantify distinct intermolecular interactions coming from different atomic pairs [44]. As shown in Fig. S5, two sharp spike-like structures were observed around the lower left of the plots, which were typical of O----H/ H----O interactions and comprised 21.9% and 25.1% of the total Hirshfeld surfaces of 1 and 2, respectively. The H----H interatomic contacts which were discovered to be another dominating interactions contributing 45.2% and 31.2% to the total surface contribution in 1 and

2.

The fingerprint plots also show the characteristic -wing‖ like structures that arise from C----H/ H----C interactions and contribute 17.7% in complex 1 and 30.2% in complex 2 of the total surface area.

Stability studies in solution

The stabilities of the complexes have been confirmed by absorption spectroscopic technique. The absorption spectra of both the complexes were performed in DMSO as well as in water at distinct time intervals of 0 h, 1 h, 12 h and 24 h. The results showed that neither the intensity nor the position of the absorption bands had significantly changed. Furthermore, even after 24 hours, no discernible changes in the absorption spectra were seen, demonstrating the stability of their geometries and coordinating spheres. (Fig. 6, S9)

ACCEPTED MANUSCRIPT / CLEAN COPY

Binding profile with therapeutic targets, ct-DNA/tRNA

In vitro interactions with complexes 1 and 2 were performed to investigate their binding mode, and their preferential selectivity towards ct-DNA and RNA as therapeutic targets. Since most of the prominent metal-based anticancer drugs in clinical use have their pharmacological target as DNA, targeting DNA provides a rationale therapeutic approach for their drug design. In the recent past, RNAs have gained paramount importance for designing efficient metal-based drug candidates due to their biological role and diverse structural motifs providing unique binding pockets for selective specific recognition. Despite this, there are fewer reports that describe the binding profile of small metal-based drug candidates such as Schiff base complexes with tRNA.

Thus, investigating the molecular mechanism of metallo-drugs targeted to tRNA biomolecule serves as a prerequisite to design efficient metal-based chemotherapeutic agents.

Absorption spectral studies

To determine the binding propensity and target selectivity of complexes 1 and 2, absorption titration experiments with ct-DNA and tRNA were performed. Titration curves of complex 1 with increasing ct-DNA/tRNA concentration in aqueous Tris-saline buffer revealed an increase in the absorption spectrum at the ligand centered π-π* (λ = 277) and while a decrease at the ACCEPTED MANUSCRIPT / CLEAN COPY changes implicate a partial intercalative binding mode of the complex 1 [45]. On the other hand, complex 2 exhibited a decrease in in the absorption intensity at π-π* (λ = 304) band and at the LMCT band (λ = 403 nm). The -hypochromic effect‖ in absorption spectra of complex 2 could be ascribed to π-π stacking interactions between the planar aromatic chromophore of Schiff base ligand scaffold and nucleobases of ct-DNA and tRNA molecules, which is compatible with an intercalative binding mode [46].

Complex 1 suggested 15.43% and 17.85% hypochromism, whereas complex 2 resulted in 19.04% and 21.41% hypochromism for ct-DNA and tRNA, respectively. The greater extent of hypochromism towards the biomolecules for 2 resulted in higher binding affinity in comparison to 1. To further quantify the extent of binding, Wolfe-Shimer equation [47] was used to calculate K b values, and the values are presented in Table 1. The values for Gibb's free energy were also evaluated for 1 and 2. Furthermore, higher values of K b upon interaction with tRNA clearly indicate the preferential selective binding of the complexes towards RNA as compared to ct-DNA.

ACCEPTED MANUSCRIPT / CLEAN COPY

Fluorescent intercalator displacement assays

To gain a better understanding of the mode of interaction, a fluorescent intercalator displacement assay with increasing aliquots of 1 and 2 using ethidium bromide was carried out. With the subsequent addition of 1 and 2 (Fig. S6), the emission intensity decreased as the bound intercalator, ethidium bromide was gradually replaced from the ct-DNA/tRNA nucleobases. The observed shifts in the emission intensity of 1 and 2 suggested an intercalative binding mode towards ct-DNA and tRNA [48]. The classical Stern-Volmer equation was applied to calculate the quenching constant Ksv values (Fig. 8) for complexes 1 and 2. The complex 2 observed higher K sv value, as shown in Table 1, which revealed its stronger quenching ability than complex 1 towards ct-DNA/tRNA biomolecule.

ACCEPTED MANUSCRIPT / CLEAN COPY

Circular dichroic studies

Circular dichroic spectral studies were carried out to analyze the conformational alterations in the secondary structure of ct-DNA and tRNA macromolecules. The CD spectrum of the B-form ct-DNA exhibited a positive peak centered at 275 nm arising from nucleobase stacking and a negative peak at 245 nm due to helicity. Intercalative binding resulted in significant perturbations of these optical signatures; however, electrostatic interactions or groove binding causes less or no disruption in these bands.

Upon addition of 1 and 2, there is an upward shift in the positive signal and a decrease in the negative signal along with a wavelength shift in the ct-DNA's CD spectrum was observed, as shown in Fig. 9. Such changes were suggestive of the intercalative mode of binding of the complexes within the DNA helix, causing unwinding and inducing conformational changes in the B-DNA structure [49].

ACCEPTED MANUSCRIPT / CLEAN COPY Upon incubation of 1 and 2, a marked increase was observed in the positive bands and a decrease in the signal intensity of the negative bands accompanied by a wavelength shift (Fig. 10). These spectral changes suggested that the complexes interact via intercalative mode resulting in conformational changes within the tRNA molecule [50].

Fig. 10. The observed shifts in CD spectra of tRNA molecule upon interaction with complexes 1 (pink) and 2 (brown). Spectrum of tRNA alone was indicated with blue line.

EPR interaction studies

Previous literature reports have shown that the copper(II) ion could attack the DNA base pairs and coordinate monodentately with the ring nitrogen atom of the guanine nucleobases at pH=7.5 [START_REF] Lippert | Coordinative bond formation between metal ions and nucleic acid bases[END_REF]. With the addition of nucleic acids, EPR signal intensity either diminutes or increases, which suggested interaction of metal complexes with the ct-DNA/tRNA. EPR spectral measurements were performed to investigate the binding of Cu(II) complexes with ct─DNA/tRNA. The low-temperature X-band EPR spectrum of 1 and 2 with ct-DNA and tRNA ACCEPTED MANUSCRIPT / CLEAN COPY at physiological pH are presented in Fig. 11. The EPR spectra of 1 and 2 in the frozen solution state revealed a single resonance at g∥ and g⊥, which were characteristic of coupled Cu(II) ion.

However, upon addition of an equivalent amount of ct-DNA/ tRNA to 1 and 2, the spectral signatures revealed a shift in EPR spectral signatures along with the hyperfine splitting having g∥ and g⊥ values corresponding to the copper complex with apical chelation by a donor atom.

The structural features of complexes 1 and 2 show that Cu(II) was coordinated through two nitrogen atoms and two oxygen atoms in the equatorial plane, with one of the oxygen atoms coordinated at the axial site. After titration with ct-DNA /tRNA, EPR spectra suggested that the Cu(II) was bound equatorially with one nitrogen from DNA/RNA base after the replacement of axially coordinated oxygen atom [START_REF] Facchin | [END_REF]. Furthermore, the signal intensity was diminished to a larger extent with the addition of tRNA, resulted in a stronger interaction of the complexes with tRNA as compared to ct-DNA. The changes observed in the EPR spectral parameters of complexes 1 and 2 in the absence and presence of ct─DNA /tRNA are shown in Table 2. The DNA cleavage mechanism plays an important role in cellular biology as it can cause DNA damage, prevent replication, and ultimately results in cell death. [53] Metal complexes could exert their cytotoxic effects by inducing DNA damage through redox activity or hydrolysis of phosphate sugar backbone. Thus, studying the cleavage mechanism provides a pathway for designing efficient metallonucleases which were capable of inducing DNA damage.

Gel electrophoresis is a frequently used technique for evaluating DNA cleavage behaviour by measuring the variations in DNA mobility when an electric field is applied to incubated complexes. When plasmid DNA is subjected to electrophoresis, the fastest travelling supercoiled form (Form I) is converted to a relaxed nicked form (Form II), implying single-strand DNA scission. Between Form I and Form II, on the other hand, the development of the linearized form of DNA (Form III) predicted more lethal double-strand breaks of the DNA helix inducing cell death.

To investigate the DNA cleaving activity of 1 and 2, a concentration-dependent pBR322 DNA agarose gel electrophoretic assay was performed under physiological conditions. With the progressive addition of complexes 1 and 2, the plasmid DNA was incubated for 1 h, and the electrophoretic bands showed that the supercoiled form (Form I) was converted to nicked circular form (Form II) along with the appearance of the linearized form (Form III) (Fig. 12 a & b). These results demonstrated that complexes 1 and 2 had the ability to induce double-strand DNA cleavage at low micromolar concentrations. Furthermore, complex 2 was able to induce ACCEPTED MANUSCRIPT / CLEAN COPY much more efficient cleavage at 5 µM as compared to complex 1 whose cleavage activity was observed at 8 µM. resulting in the formation of copper-oxo species [54,55].

To investigate the complexes' preferential groove binding mode, cleavage activity was measured with the addition of minor groove binding agent, DAPI (4′,6-diamidino-2-phenylindole). The electrophoretic bands (Fig. 13 a & b, lane 8) demonstrated that the addition of DAPI resulted in considerable cleavage inhibition activity resulted in the minor groove preference of 1 and 2.

These findings are of considerable interest as most of the oxidative cleaving reagents interact preferentially through minor grooves rather than major grooves [56].

RNA Cleavage studies

ACCEPTED MANUSCRIPT / CLEAN COPY

The RNA molecule was involved in various cellular processes and performed regulatory functions, thus could be exploited as a viable therapeutic target. The metal-based complexes were found to promote specific degradation of the RNA molecule through hydrolysis of the phosphate backbone or as a sequence-recognizing moieties. Thus, investigating the RNA cleavage pathway in the presence of metal complexes could be envisaged for controlling RNA structural features and developing metallonucleases that could selectively hydrolyze RNA at the target site.

The RNA cleavage activity was monitored at a fixed concentration of tRNA which was incubated with an increasing concentration of the complexes 1 and 2 at room temperature for a period of 3 h and 24 h. After 3 hours of incubation, the complexes caused a partial cleavage of the tRNA molecule and showed moderately decreased bands, according to the gel electrophoretic pattern (Fig. 14 a & b). However, with concentration range from 6.25 µM to 17.50 µM, the cleavage activity was increased and found to be higher in complex 2 than 1. The cleavage activity was also monitored after 24 h of incubation which revealed the occurrence of largely diminished electrophoretic bands suggesting complete degradation (90-95%) of the tRNA (Fig. 14 c &d). Furthermore, with increasing concentration, the complexes demonstrated higher activity implicated by the diminishing band intensity, which was higher in complex 2. These findings showed that RNA cleavage was time-and concentration-dependently promoted by complexes 1 and 2 [57].

Lipophilicity

Lipophilicity is used to determine various biological parameters of the complexes, such as cellular uptake, metabolism, bioavailability, and lipid solubility. The study has gained importance as it is ideal for predetermining the ADMET properties of drug candidate correlate structure and lipophilicity profiles, thus, rendering them a safe throughput in clinical trials. The lipophilicity can be calculated in terms of Log P, demonstrating the membrane permeability of the complexes. The SwissADME software calculated Log P values of complexes 1 and 2 as 1.28 and 2.82, respectively. The greater lipophilic nature of complex 2 than complex 1 was reflected by their LogP value, followed by their enhanced biological activity [58].

Cytotoxic profile

The cytotoxicity of both the complexes was assessed using MTT assay against three cancer cell lines viz., MDA-MB-231 (triple negative breast cancer cell line), HCT 116 (human colon cancer cell line) and A549 (lung cancer cell) (Fig. 15). When 1 and 2 were treated at a dosage of 10-200 µg/ml, it was found that the complexes effectively inhibited the growth of these malignant strains in a dose-dependent manner. The complex 2 observed greater inhibitory response as quantified by the calculated IC 50 values (Table 3). This is caused due to its unsaturated and planar geometric preference that facilitates the insertion of 2 into the plasma membrane. The complex displayed a stronger inhibitory response on the tested MDA-MB-231 cancer cell line as a result of the positive association between DNA/RNA interaction and cytotoxicity. Furthermore, ROS species produced by the complexes may be responsible for inducing apoptosis in cancer cells by producing H 2 O 2 levels above the toxic threshold. In vitro tumour cell growth is directly and dose-dependently inhibited when anticancer drugs are combined with ROS species. The morphological alterations detected in the cell lines exposed to complexes for 48 h are presented in Fig. 16 indicating that with the exposure of complexes to the cell lines, complexes 1 and 2

ACCEPTED MANUSCRIPT / CLEAN COPY morphological changes indicated that these complexes are capable of inducing cell death, probably by apoptosis cells. However, changes in complex 2 were more pronounced as compared to complex 1, making it as a potential cytotoxic agent. ACCEPTED MANUSCRIPT / CLEAN COPY Fig. 16 The morphological alterations in different cancer cells upon treatment with complexes at their IC 50 values (µg/ml)

Molecular docking studies

Molecular docking studies with DNA and RNA biomolecule were conducted to further confirm the binding mode as well as selective recognition sites for 1 and 2, were carried out [59]. The Molecular docking studies were also performed with DNA and RNA biomolecules which revealed higher negative binding energy values of the complexes with tRNA, validating the results obtained from spectroscopic studies.

Complex [C 13 H

 13 15 CuN 3 O 5 S], 1: Yield: 85%, M.P: 235˚C, Anal. (%) Calc. for C 13 H 15 CuN 3 O 5 S: C (41.73), H (4.25), N (10.43), Found: C (41.74), H (4.24), N (10.45). FT-IR (KBr cm -1 ): 1622 ν(HC=N), 1243, 1139 ν(S=O), 1474 ν(C=N) imidazole, 543 ν(M-O), 486 ν(M-N). UV-vis (λ max , nm) in DMSO: 263 (π -π*), 347 (n-π*), 678 (d-d). ESI-MS (m/z): 387 [C 13 H 14 CuN 3 O 5 S] + . Complex [C 16 H 15 CuN 3 O 4 S], 2: Yield: 87%, M.P: 248˚C, Anal. (%) Calc. for C 16 H 15 CuN 3 O 4 S: C (46.88), H (3.93), N (10. .43), Found: C (41.74), H (4.24), N (10.45). FT-IR (KBr cm -1 ): 1618 ν(HC=N), 1236, 1136 ν(S=O), 1434 ν(C=N)imidazole, 588 ν(M-O), 499 ν(M-N). UV-vis (λ max , nm) in DMSO: 266 (π -π*), 354 (n-π*), 680 (d-d). ESI-MS (m/z): 403 [C 16 H 14 CuN 3 OS-3H]+.

complex 1

 1 showed a molecular ion peak at m/z = 387 corresponding to the [C 13 H 14 CuN 3 O 5 S] + and an isotopic peak at m/z = 389. Another peak was observed due to the removal of imidazole ligand moiety at m/z = 319 attributed to the molecular fragment [C 10 H 10 CuNO 5 S] + along with an isotopic peak at m/z = 320. Complex 2 revealed a peak at m/z = 403 corresponding to the fragment [C 16 H 14 CuN 3 O 4 S-3H + ] and other peak at m/z = 340 due to removal of imidazole ligand with [C 13 H 10 CuNO 4 S+2H + ]. Complexes 1 and 2 also showed peaks at m/z = 257 and 278, which are assignable to the ligand moieties with molecular fragments [C 10 H 12 NO 5 S-H] + and [C 13 H 12 NO 4 S] + . Thus, the mass spectra of complexes 1 and 2 correlated well with their predicted molecular structures.

  possesses a longer Cu-O bond distance (2.506 Å) than Cu-O bonds (1.907-2.001Å) in the basal ACCEPTED MANUSCRIPT / CLEAN COPY plane, demonstrating the presence of Jahn-Teller distortion. The elongation results from the unsymmetrical occupation of the electrons present in the d-orbital of the Cu(II) atom, where two electrons are found in the d z 2 while one electron is accommodated in the d deviation from the equatorial plane in complex 1 was also ascertained by measuring the trans angles, i.e., N(1)-Cu(1)-O(19) and N(11)-Cu(1)-O(31), which were found to be 88.33º and 94.02º, respectively. The trigonality index (τ) value was calculated according to Addison-Reedjik geometric parameter i.e., τ = (β-α)/60 where α and β represents the two largest angles around the copper center. The value of trigonality index for ideal square pyramidal geometry is zero and for perfect trigonal bipyramidal τ = 1. In our case, at the Cu(II) site it was found to be τ = 0.069, close to the index expected for a distorted square pyramidal geometry. The two Cu(II) metal centers are 5.10 Å distance apart from one another. Significant intermolecular hydrogen and C---C bonding interactions were also present along the c axis in a one-dimensional polymeric array of the crystal (Fig. 2). Complex 1 showed extensive hydrogen bonding between the hydrogen atoms of imidazole and oxygen atom of the Schiff base ligand and C---C bonding interactions between the C5 carbon atoms of imidazole rings arranged in an antiparallel fashion generating a one dimensional (1D) polymeric array along the c axis.

Fig. 1 .

 1 Fig. 1. The crystal structure of complex 1 drawn at 50% probability level.

Fig. 2 .

 2 Fig. 2. 1D polymeric structure in complex 1 showing intermolecular (a) hydrogen and (b) carbon-carbon bonding interactions along the c axis.The crystallographic details of complex 2 revealed that the complex crystallizes in a monoclinic

Fig. 3 .

 3 Fig. 3. ORTEP view for the crystal structure of complex 2 drawn at 50% probability level.

Fig. 4 .

 4 Fig. 4. 2D ladder-like polymeric array in complex 2 revealing intermolecular hydrogen bonding interactions along the c axis.

  and the HOMO and LUMO energy gap was associated with the kinetic stability and chemical reactivity. The analysis of FMOs of 1 and 2 demonstrated that the LUMOs mainly outspread over the amine nitrogen and the aromatic ring of the Schiff base ligand while the HOMOs were dispersed around the metal coordinated sulfonate group of the Schiff base and the imidazole ACCEPTED MANUSCRIPT / CLEAN COPY ligand. The LUMOs are responsible for the electron-donating ability of the coordinated ligand in the unoccupied metal-centered orbitals, whereas the HOMOs show the electron-accepting ability

Fig. 5 .

 5 Fig. 5. Frontier molecular orbitals contour plot of 1 and 2.

Fig. 6 .

 6 Fig. 6. Stability studies of complexes 1 (a) and 2 (b) in DMSO at different time intervals (0h, 1h, 12h, 24 h), [Complex 1] = [Complex 2] = 10 -3 M

Fig. 7 .

 7 Fig. 7. Changes in absorption intensity of complexes 1 and 2 with progressive addition of DNA (a and b) and tRNA (c and d), respectively.

Fig. 8 .

 8 Fig. 8. K sv plots of complexes 1 and 2 towards binding with (a) ct-DNA and (b) tRNA..

Fig. 9 .

 9 Fig. 9. CD spectra of ct-DNA alone (blue) and in the presence of complexes 1 (pink) and 2 (brown). Due to ellipticity, the CD spectrum of tRNA revealed two positive bands at 270 and 221 nm,

Fig. 11 .

 11 Fig. 11. EPR spectra of 1 and 2 (black) show shifts in the signal amplitude on interaction with ct-DNA (red) and tRNA (green).

Fig. 12 .

 12 Fig. 12. Complexes 1 and 2 showed pBR322 DNA cleavage (0.1 μg/μL) at different concentrations at 37 °C with an incubation time of 1 h: (a) lane 1, DNA; lane 2, DNA + 1 (1μM); lane 3, DNA + 1 (2μM); lane 4, DNA + 1 (3μM); lane 5, DNA + 1 (4μM); lane 6, DNA + 1 (5μM); lane 7, DNA + 1 (6μM); lane 8, DNA + 1 (7μM); lane 9, DNA + 1 (8μM); lane 10, DNA + 1 (9μM); lane 11, DNA + 1 (10μM); lane 12, DNA + 1 (11μM); (b) lane 1, DNA; lane 2, DNA + 2 (1μM); lane 3, DNA + 2 (2μM); lane 4, DNA + 2 (3μM); lane 5, DNA + 2 (4μM); lane 6, DNA + 2 (5μM); lane 7, DNA + 2 (6μM); lane 8, DNA + 2 (7μM); lane 9, DNA + 2 (8μM); lane 10, DNA + 2 (9μM); lane 11, DNA + 2 (10μM); lane 12, DNA + 2 (11μM).Earlier reports observed that the copper(II) complexes had the ablilty to cleave the plasmid DNA

Fig. 13 .

 13 Fig. 13. cleavage behavior of pBR322 DNA (0.1 μg/μL) in the prescence of 1 and 2 at 37 °C with an incubation time of 1 h with the addition of different additives: (a) lane 1, DNA; lane 2, DNA + 1; lane 3, H 2 O 2 (0.25 μM) + DNA + 1; lane 4, sodium azide (0.4 μM) + DNA + 1; lane 5, SOD (0.4 μM) + DNA + 1; lane 6, Ethyl alcohol (0.4 μM) + DNA + 1; lane 7, DMSO (0.4 μM) + DNA + 1; lane 8, DAPI (0.16 μM) + DNA + 1; (b) lane 1, DNA; lane 2, DNA + 2; lane 3, H 2 O 2 (0.25 μM) + DNA + 2; lane 4, sodium azide (0.4 μM) + DNA + 2; lane 5, SOD (0.4 μM) + DNA + 2; lane 6, Ethyl alcohol (0.4 μM) + DNA + 2; lane 7, DMSO (0.4 μM) + DNA + 2; lane 8, DAPI (0.16 μM) + DNA + 2.

Fig. 14 .

 14 Fig. 14. The gel electrophoresis pattern of tRNA molecule with the addition of increasing concentration of complexes 1 and 2 (6.25-17.50 µM) for 3h (a and b) and for 24h (c and d).

Fig. 15

 15 Fig. 15 Cell viability of treated cancer cells at different complexes concentrations under dark and

docked pose model of 1 and 2

 2 with DNA (PDB ID: 1BNA) delineated that the complexes preferentially target the minor groove GC region in proximity to G-24, C-23, G-22, G-4, G-4, A-5 and A-6 nucleobase pairs (Fig. S7 & S8a). The Schiff base ligand in the complexes was found to intercalate within the DNA nucleobases, while the imidazole moiety interacted through hydrogen bonding interactions stabilizing the overall molecular interactions. The binding energies were found to be -248 KJmol -1 and -270 KJ mol -1 for complexes 1 and 2, respectively, suggesting the higher binding propensity of complex [C 16 H 15 CuN 3 O 4 S] with DNA duplex. ACCEPTED MANUSCRIPT / CLEAN COPY 1. Furthermore, the complexes show higher binding propensity towards tRNA as compared to ct-DNA implicating their selective preferential binding towards RNA. The gel electrophoretic experiments of complexes 1 and 2 were carried out with plasmid pBR322 DNA which demonstrated that the complexes were able to cleave the plasmid DNA in a dose dependent manner at low micromolar concentrations. The mechanistic insight of DNA cleavage in presence of ROS radical scavengers showed that DNA scission was facilitated suggesting the involvement of singlet oxygen ( 1 O 2 ) and superoxide anions (O 2 •-). The cleavage studies in presence of H 2 O 2 implicated a lethal double-strand DNA scission via the oxidative pathway involving the formation of copper-oxo species mediated through Fenton's reaction mechanism. The tRNA cleavage experiments of the complexes 1 and 2 were performed which revealed efficient cleavage activity of the complexes as a function of both increasing concentration and time.

Table 1 .

 1 Binding values of complexes 1 and 2 on the interaction with ct-DNA and tRNA.

Complexes

K b (M

-1 ) ΔG (KJ/mol) % Hypo chromism K sv (10 5 M -1 )

		DNA	RNA	DNA	RNA	DNA	RNA	DNA	RNA
	1	4.37(±0.04) x10 5 5.62(±0.02) x10 5 -32.15 -32.77 15.43 17.85	1.07	1.37
	2	5.82(±0.05) x10 5 6.34(±0.03) x10 5 -32.87 -33.07 19.04 21.41	1.15	2.21

Table 2 .

 2 EPR spectral parameters of complexes 1 and 2.

	Complexes	Complex alone	Complex + ct-DNA	Complex + tRNA
	1	g || =2.34, g⊥=1.94	g || =2.47, g⊥=1.88	g || =2.32, g⊥=1.89
	2	g || =2.39, g⊥=1.93	g || =2.45, g⊥=1.90	g || =2.31, g⊥=1.88
	3.9.5 Cleavage studies pBR322 DNA

Table 3 Complexes IC 50 values (µg/ml)

 3

		MDA-MB-231	HCT-116	A549
	1	118.5	81.97	77.61
	2	22.96	51.59	53.53

Acknowledgements

The SAIF Panjab University, Chandigarh is gratefully acknowledged for providing elemental analysis and ESI-MS Facilities. We are also thankful to the Department of Chemistry, Aligarh ACCEPTED MANUSCRIPT / CLEAN COPY

The Docking profile of the complexes 1 and 2 with tRNA (PDB ID: 6tRNA) demonstrated that these complexes fit snugly into the RNA intercalation site i.e., inserted into the active pocket of the anticodon arm in proximity to C-28, C-27, A-29, G-26, C-25, G-30, C-40, G-42, G-43 and U-41 (Figure S7 & S8b). In addition, van der Waals and hydrophobic interactions also play a major role in the stronger binding of the complexes with tRNA. The binding energies for complexes 1 and 2 were found to be -265 KJmol -1 and -299 KJmol -1 , respectively, which validates the higher binding propensity of complex 2. Furthermore, the higher negative binding energy values of the complexes with tRNA as compared to ct-DNA demonstrated that they were better targeted towards tRNA therapeutic biomolecule.

Conclusions

New water soluble polymeric (μ-oxo) bridged Cu(II) complexes derived from taurine Schiff base and appended with imidazole ligand were synthesized. The structures of complexes were established by various spectroscopic and single X-ray crystallographic studies. The complex 1 crystallizes in the trigonal crystal system with the R-3 space group and complex 2 crystallizes in a monoclinic crystal system with the P 21/n space group. The crystal structure of the complexes was further stabilized by the intermolecular C----C and H----H bonding interactions resulting in a two dimensional supramolecular framework.

The DFT studies were carried out to determine the frontier molecular orbitals and other geometrical parameters of the complexes and to correlate their energy and stability. The results of DFT studies demonstrated higher stability of complex 2 as compared to 1.

The binding profile of complexes

Conflict of Interest

There is no conflict of interest to declare.


Javier Antonio Gonzalez

Laurent Dupuy

Peter Råback

Marc Fivel

Michel Perez

Jonathan Amodeo
email: jonathan.amodeo@cnrs.fr

El-Numodis: a new tool to model dislocation and surface interactions

While surfaces are known to have a limited impact on the mechanical properties of crystalline materials at the macroscopic scale, they play a key role at small-scale behaving alternatively as sources or sinks of various plastic deformation processes. In this study, we present a new tool called El-Numodis that relies on the superposition method to couple the discrete dislocation dynamics code Numodis to Elmer, an open-source finite-element-modeling tool. After few years of development, El-Numodis allows now for the simulation of small-scale object deformation and mechanical properties based on a large set of surface-related processes including stress-free boundaries, mirrored dislocations and a Monte-Carlo based dislocation nucleation mechanism. Here we present the main features of the code as well as numerical test-cases and benchmarks going from classical boundary value problems to tensile tests on model thin film.

Indeed, nanotech is now an important segment of the modern material industry. The wide range of applications of nanodevices is mainly due to the advantages brought by the radical change of their properties (mechanical, surface, optical, chemical, etc.) induced by their size, compared to bulk materials. As examples, nanowires are employed for digital data storage due to their superparamagnetic capability that relies on a very fast response to external fields with almost zero remanence [START_REF] Hui Lu | Magnetic nanoparticles: Synthesis, protection, functionalization, and application[END_REF][START_REF] Zhang | Synthesis and assembly of magnetic nanoparticles for information and energy storage applications[END_REF], nanoparticles are also currently used to improve the performance of lubricants [START_REF] Akbulut | Nanoparticle-Based Lubrication Systems[END_REF][START_REF] Dai | Roles of nanoparticles in oil lubrication[END_REF] or as compounds to build implants due to their high-strength and wear resistance [START_REF] Li | Applications of nanotechnology in hip implants[END_REF]. Overall, nanocrystals and their outstanding mechanical properties (smaller is stronger ) are now widely used to improve bulk materials [START_REF] Taneike | Creep-strengthening of steel at high temperatures using nano-sized carbonitride dispersions[END_REF][START_REF] Casati | Metal matrix composites reinforced by nano-particles-a review[END_REF] and their strength vs. size dependence mostly arises from the increase of the surface over volume ratio when scaling down the sample size [START_REF] Greer | Plasticity in small-sized metallic systems: intrinsic versus extrinsic size effect[END_REF][START_REF] Mordehai | Nucleation-controlled plasticity of metallic nanowires and nanoparticles[END_REF]. While the flow of bulk materials is known to be governed by a dislocation multiplication process from an existing defective microstructure, the plastic deformation of nano-objects is controlled by a Surface Dislocation Nucleation (SDN) mechanism that requires a much larger stress than in the classical dislocation multiplication case (stress increased by a factor of 10 to 1000) [START_REF] Andrew | Heterogeneous dislocation nucleation from surfaces and interfaces as governing plasticity mechanism in nanoscale metals[END_REF][START_REF] C R Weinberger | Atomistic simulations and continuum modeling of dislocation nucleation and strength in gold nanowires[END_REF][START_REF] Amodeo | Mechanical properties and dislocation nucleation in nanocrystals with blunt edges[END_REF][START_REF] Issa | From dislocation nucleation to dislocation multiplication in ceramic nanoparticle[END_REF]. In fact, nanocrystals are known to be dislocation-scarce (or free) due to both the soft fabrication routes from which they are derived (e.g., crystal growth, dewetting, lithography) as well as to surface-induced image forces [START_REF] Hirth | Theory of dislocations[END_REF][START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF] that pump the dislocation density out, this latter being intrinsically concomitant to the aforementioned SDN process.

Several experimental and numerical methods currently exist to try to understand the mechanical behaviour of nano-objects. For example, nanocompression in the Scanning or Transmission Electron Microscope (SEM or TEM) are the most used experimental techniques in the field [START_REF] Minor | Quantitative in situ nanoindentation in an electron microscope[END_REF][START_REF] Östlund | Brittle-to-ductile transition in uniaxial compression of silicon pillars at room temperature[END_REF][START_REF] Legros | Quantitative In Situ Mechanical Testing in Electron Microscopes[END_REF]. Experiments are reported for sample sizes ranging from several micrometers down to few tens nanometers as well as for low deformation rate (10 -4 to ∼1 s -1 ). While the use of TEM allows for microstructure, defect and surface characterizations at the nanoscale, it is still a complex and expensive method and several issues including sample misalignment, contamination or oxidation are commonly reported [START_REF] Williams | The Transmission Electron Microscope[END_REF][START_REF] Zaefferer | A critical review of orientation microscopy in SEM and TEM[END_REF][START_REF] Gerberich | Review Article: Case studies in future trends of computational and experimental nanomechanics[END_REF]. Thus, computational methods are often used to support the interpretation of nanomechanics experiments. On one hand, MD is the most used numerical technique to simulate the mechanical properties of nano-objects as relying on the description of atomic-scale processes [START_REF] Mordehai | Nucleation-controlled plasticity of metallic nanowires and nanoparticles[END_REF][START_REF] Amodeo | Modeling the mechanical properties of nanoparticles: a review[END_REF]. It is based on interatomic potentials used to compute atomic forces and integrate the dynamics of molecular systems e.g., in a sample under load when applied to the field of nanomechanics. However, MD sample size is usually limited to few tens of nanometers while nano-objects can be of few hundreds and, furthermore, MD is performed at particularly high strain-rate (∼10 8 s -1 ) due to prohibitive computational costs. On the other hand, while DDD allows for larger sample size and strain rates both being closer to experimental conditions, it also has some drawbacks. Indeed, DDD allows to model dislocations on the basis of constitutive equations rather than atomic interactions what makes the model less expensive in terms of cpu costs but has the disadvantage of being less accurate. Furthermore, DDD codes are generally developed for bulk applications and often miss nanoscale features. The major lack relies on the elastic theory the DDD is based on that usually assumes an infinite continuum. This peculiarity raises the issue of finite Boundary Conditions (BCs) as DDD was originally developed to tackle configurations close to bulk conditions i.e., using Periodic Boundary Conditions (PBCs) or free-BCs, without explicitly considering realistic surfaces. In addition, former DDD codes (e.g., microMegas [START_REF] Devincre | The modelling of dislocation dynamics: elastic behaviour versus core properties[END_REF][START_REF] Devincre | Modeling crystal plasticity with dislocation dynamics simulations: the 'microMegas' code[END_REF], Paradis [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF], Tridis [START_REF] Fivel | 3D simulation of a nanoindentation test at a mesoscopic scale[END_REF][START_REF] Marc | Discrete dislocation dynamics: an important recent break-through in the modelling of dislocation collective behaviour[END_REF]) have not been originally developed to describe dislocation nucleation processes. While this latter is not crucial to simulate bulk mechanics, it is imperative to model nanoscale mechanical properties due to the aforementioned SDN process. In this context, the development of quantitative numerical tools to investigate the mechanical properties of small-scale objects including the effect of surfaces, sample size, strain rate as well as the ability to account for most of nanoscale elementary deformation processes emerges as a crucial step up to improve our understanding of micro and nanomechanics [START_REF] Gerberich | Review Article: Case studies in future trends of computational and experimental nanomechanics[END_REF].

Boundary Value Problems (BVPs) related to surfaces or interfaces can be addressed coupling the DDD with the Finite-Element Method (FEM) using various kinds of approaches such as the SuperPosition Method (SPM) that relies on a linear stress correction originally proposed by Van der Giessen and Needleman [START_REF] Van Der Giessen | Discrete dislocation plasticity: A simple planar model[END_REF], the Discrete Continuous Method (DCM) proposed by Devincre and collaborators [START_REF] Lemarchand | Multiscale modelling of plastic deformation[END_REF][START_REF] Vattré | Modelling crystal plasticity by 3D dislocation dynamics and the finite element method: The Discrete-Continuous Model revisited[END_REF][START_REF] Jamond | Consistent formulation for the Discrete-Continuous Model: Improving complex dislocation dynamics simulations[END_REF] based on the eigenstrain formalism of Mura [START_REF] Mura | Micromechanics of Defects in Solids[END_REF] or the Fast Fourier Transform (FFT) dislocation approach which is increasingly used by the community [START_REF] Bertin | A FFT-based formulation for efficient mechanical fields computation in isotropic and anisotropic periodic discrete dislocation dynamics[END_REF][START_REF] Gao | Influence of misfit stresses on dislocation glide in single crystal superalloys: A three-dimensional discrete dislocation dynamics study[END_REF]. As an example of application, a DCM-FFT approach was recently used by Kohnert and colleagues to quantify the effect of surfaces and TEM lamella thinning on the dislocation density [START_REF] Aaron A Kohnert | On the use of transmission electron microscopy to quantify dislocation densities in bulk metals[END_REF].

Here we present our approach called El-Numodis which is based on the SPM and integrate specific features adapted to nanomechanical simulations. El-Numodis relies on the coupling of the DDD code Numodis [START_REF] Drouet | Dislocation dynamics simulations of interactions between gliding dislocations and radiation induced prismatic loops in zirconium[END_REF][START_REF] Shi | Interaction of < 100 > dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF] and the open-source FEM software Elmer [START_REF] Malinen | Elmer finite element solver for multiphysics and multiscale problems[END_REF]. It benefits of a particularly accurate and parallel nodal DDD framework that integrates top of the art features such as various implementations of the elastic theory, the singular and non-singular theories for dislocations [START_REF] Hirth | Theory of dislocations[END_REF][START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] as well as additional ingredients to better model deformation tests at small-scales accounting for surface effects as e.g., the Weygand's approach for dislocation vs. surface interactions [START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF] and a Monte-Carlo routine for homogeneous and heterogeneous dislocation nucleation, as inferred from the harmonic Transition State Theory (TST) [START_REF] Zhu | Temperature and strain-rate dependence of surface dislocation nucleation[END_REF]. In the following, we start with a brief reminder about the DDD and FEM techniques in the context of the parent codes Numodis and Elmer, then, the main aspects of the coupling are introduced (SPM basics, coupling and interfacing procedures, dislocation nucleation and interaction with surfaces, etc.). Finally, an extended last section composed of various validation tests and applications complete the study.

Numerical methods

Basics on DDD and FEM: methodologies and parent codes

The DDD is a mesoscopic approach developed to investigate the collective behaviour of dislocations. This method was originally developed in the late 80' to study the evolution of a dislocation population under load in metals [START_REF] Lepinoux | The dynamic organization of dislocation structures: a simulation[END_REF]. Since, it was further used to investigate various issues including strain hardening (see e.g., [START_REF] Vv Bulatov | Dislocation multi-junctions and strain hardening[END_REF][START_REF] Madec | From dislocation junctions to forest hardening[END_REF]), irradiation defects [START_REF] Shi | Interaction of < 100 > dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF][START_REF] Diaz De La Rubia | Multiscale modelling of plastic flow localization in irradiated materials[END_REF][START_REF] Arsenlis | A dislocation dynamics study of the transition from homogeneous to heterogeneous deformation in irradiated body-centered cubic iron[END_REF], fatigue [START_REF] Déprés | Low-strain fatigue in AISI 316L steel surface grains: a three-dimensional discrete dislocation dynamics modelling of the early cycles I. Dislocation microstructures and mechanical behaviour[END_REF][START_REF] Gagel | Formation of extended prismatic dislocation structures under indentation[END_REF] or nanoindentation [START_REF] Hyung | Multiscale modelling of nanoindentation[END_REF] in metals. Few applications in oxides and minerals can also be found [START_REF] Portelette | Athermal dislocation strengthening in uo 2[END_REF][START_REF] Amodeo | Multiscale modelling of MgO plasticity[END_REF][START_REF] Amodeo | Dislocation reactions, plastic anisotropy and forest strengthening in MgO at high temperature[END_REF][START_REF] Reali | Modeling viscosity of (Mg,Fe)O at lowermost mantle conditions[END_REF]. General details about the DDD method can be found in Refs. [START_REF] Devincre | Modeling crystal plasticity with dislocation dynamics simulations: the 'microMegas' code[END_REF][START_REF] L P Kubin | Dislocations, Mesoscale Simulations and Plastic Flow[END_REF][START_REF] Po | Recent Progress in Discrete Dislocation Dynamics and Its Applications to Micro Plasticity[END_REF]. Numodis is a versatile DDD code developed by the french Atomic Energy Center (CEA) used here in the DDD/FEM coupling framework of El-Numodis. Up to now, Numodis has shown to be particularly suited to investigate the influence of radiation-induced defects on the plastic behavior of metals such as iron [START_REF] Shi | Interaction of < 100 > dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF][START_REF] Li | Diffusion and interaction of prismatic dislocation loops simulated by stochastic discrete dislocation dynamics[END_REF] or zirconium alloys [START_REF] Drouet | Dislocation dynamics simulations of interactions between gliding dislocations and radiation induced prismatic loops in zirconium[END_REF] with one-to-one cross-validations of the DDD outcomes against MD and experiments. From a technical point of view, Numodis is a C++ nodal code in which dislocations are described by nodes interconnected into segments characterized by their Burgers vector and glide plane. It can be used either serial or parallel using the openMP protocol. In Numodis, a remeshing algorithm ensures that each segment length respects an admissible size range chosen in accordance with the characteristic size of the investigated phenomenon. The force acting on the nodes is computed within the singular [START_REF] Hirth | Theory of dislocations[END_REF] or nonsingular [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] dislocation theoretical framework using the analytical formulation derived by Arsenlis et al. [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF]. Additional dislocation core forces can be accounted using core energy terms [START_REF] Geslin | Thermal fluctuations of dislocations reveal the interplay between their core energy and long-range elasticity[END_REF][START_REF] Hu | Atomistic dislocation core energies and calibration of non-singular discrete dislocation dynamics[END_REF][START_REF] Bertin | Core energies of dislocations in bcc metals[END_REF]. The velocity of each node is computed using a classical variational approach [START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF][START_REF] Cai | A non-singular continuum theory of dislocations[END_REF] assuming an overdamped motion and various kind of dedicated mobility laws including e.g., viscous or thermally-activated glide assigned on each segment character and slip system. Dislocation contact reactions (junction, annihilation, crossedstates, etc.) are computed using the elastic theory allowing to model the dislocation microstructure evolution and strain hardening. The methodology used in El-Numodis is inspired from the seminal work of Bulatov and collaborators [START_REF] Arsenlis | Enabling strain hardening simulations with dislocation dynamics[END_REF][START_REF] Bulatov | Computer Simulations of Dislocations[END_REF] and can be resumed as follow. A collision detection algorithm is used at each timestep to predict incoming dislocation contact reactions with other dislocations or microstructural defects such as grain boundaries or precipitates. When a collision is detected, a new node is generated at the contact point and is kept fixed during the current DDD time step allowing for the rest of the dislocation microstructure to relax. If the collision is confirmed, the node evolves at step n + 1 depending on the situation while minimizing the energy using a split node algorithm. For example, for a contact between two dislocations, the node can split into two subsequent nodes forming a junction segment. At each time step, the code checks whether or not a node has to be split, a splittable node being defined as connected at least to three other nodes and must not be arbitrary pinned (e.g., Frank-Read source).

Used in the El-Numodis context as an elastic solver, Elmer is a multiphysical FEM code written in modern Fortran that includes a large set of continuum-based physical models [START_REF] Malinen | Elmer finite element solver for multiphysics and multiscale problems[END_REF]. Up to now, Elmer was used in various fields of application to model e.g., crystal growth [START_REF] V Savolainen | Simulation of large-scale silicon melt flow in magnetic czochralski growth[END_REF], blood flow in elastic arteries [START_REF] Järvinen | A method for partitioned fluid-structure interaction computation of flow in arteries[END_REF], computational glaciology [START_REF] Gagliardini | Capabilities and performance of elmer/ice, a newgeneration ice sheet model[END_REF] or electrical machines [START_REF] Keränen | Efficient parallel 3-d computation of electrical machines with elmer[END_REF]. Elmer uses high-level abstraction when treating individual equations for solving multiphysical problems and benefits from a modular structure and generic strategies that are useful when coupling it with other codes as e.g., with OpenFOAM [START_REF] Vencels | Eof-library: Open-source elmer fem and openfoam coupler for electromagnetics and fluid dynamics[END_REF]. In the following, Elmer is coupled to Numodis using the SPM method [START_REF] Van Der Giessen | Discrete dislocation plasticity: A simple planar model[END_REF][START_REF] Fivel | Implementing image stresses in a 3D dislocation simulation[END_REF]. The Elmer physical model used for this coupling is the elasticity equation that can be solved using various types of 2D (triangular, quadrilateral) and 3D (tetrahedron, hexahedron, prism or wedge) elements. Finally, the linear system is solved using number of different direct (Umfpack, MUMPS and Pardiso packages) or iterative (conjugate gradient, basic preconditioning, Krylov subsapce methods, etc.) techniques. Also, Elmer is interfaced with Hypre that provides an additional set of iterative solvers and preconditioners [START_REF] Robert | hypre: A library of high performance preconditioners[END_REF]. Finally, efficient octree-based interpolation methods that can be performed on nodes or on integration points are available for mapping results between computational meshes.

The superposition method and El-Numodis

Introduction to the superposition method

The SPM method was first introduced by Needleman and collaborators to solve BVP problems [START_REF] Van Der Giessen | Discrete dislocation plasticity: A simple planar model[END_REF]. It is based on the correction of the dislocation self-stress σ and displacement fields ũ originally computed by a DDD code (assuming an infinite medium) at physical boundaries S using a FEM elastic solver. Assuming a finite-size volume, the SPM aims at imposing an applied force or displacement (F app and u app , respectively) via the BCs that is corrected by an homologous contribution (σ or ũ) induced by a dislocation population as computed by the DDD using an infinite medium hypothesis (see Figure 1). For that purpose, the dislocation stress field computed at the boundary σ(S) is converted into forces F (S) using an appropriate conversion scheme when accounting for Neumann applied force BCs whereas applied displacement (Dirichlet) BCs ũ(S) are computed using the Barnett approach [START_REF] Barnett | The displacement field of a triangular dislocation loop[END_REF]. After solving the BVP, the field corrections ( σ or u) is added to the original internal fields (σ or ũ) computed by the DDD what leads to the total stress σ or displacement fields u (Equation 1). The FEM correction is computed at the dislocation Gauss points using a classical interpolation method and is further added to the dislocation stress or displacement field contribution. Thus, SPM allows for the mesoscale modeling of a finite-size domain including surfaces and interfaces (in the contrary to the self-standing DDD). This method that allows for arbitrary BCs was already used for various modeling applications such as nanoindendation [START_REF] Fivel | 3D simulation of a nanoindentation test at a mesoscopic scale[END_REF][START_REF] Hyung | Simulations multi-échelles de l'indentation de métaux CFC: De l'atome au milieu continu[END_REF], thin film [START_REF] Hartmaier | Image stresses in a free-standing thin film[END_REF] or micropillar compression [START_REF] Zhou | Discrete dislocation dynamics simulations of plasticity at small scales[END_REF] simulations.

σ = σ + σ u = ũ + u (1)

El-Numodis operation

In the coupling approach, the FEM code Elmer drives the simulation and refers to the DDD as an external library. Number and format conversion drivers have been implemented for direct data transfer between the two codes while Elmer external solvers were upgraded here to manage the newly coupled BCs. In a nutshell, the coupling consists in three external routines and several Numodis functions being called as shown in Figure 2. The El-Numodis workflow can be described as follows:

(i) Loading deformation conditions. El-Numodis requires a 3D geometrical mesh as well as a parameter file that are loaded at first. Among others, the parameters file contains the definition of the faces where the BCs (Dirichlet or Neumann) are applied, the total number of time steps, the output saving frequency, material elastic properties as well as the elastic solver parameters. At the end of this first step, Numodis is called using a first driver refered as El-Numodis export bnodes.

(ii) Calculation of σ(S) and ũ(S). Specific DDD inputs as e.g., the dislocation density distribution, are loaded by Numodis. Elmer associates parts of or all the surface nodes to specific BCs i.e., Neumann (including traction-free nodes) or Dirichlet BCs. Then, the DDD computes the displacement and stress fields associated with the current dislocation density σ(S) and ũ(S) at the specified BC nodes.

(iii) Fields regularisation. A second driver called El-Numodis import converts σ(S) into F (S) at Neumann boundary nodes and F (S)=F app -F (S) as well as u(S)=u appũ(S) are computed. Then, the assembly is performed accounting for a feedback loop to either adapt the applied stress or displacement. Finally, the linear elastic solver of Elmer computes σ everywhere inside the simulation volume.

(iv) Fields superposition. Finally, the last driver El-Numodis export stress is called. It first lists dislocation Gauss points and then interpolates σ. The stress superposition as described by Equation 1 is performed leading to the effective stress state σ.

Then, a new DDD step is performed including force computation, possible nucleation event or dislocation displacement. The whole process is repeated up to the total amount of simulation time steps.

Field conversion and interpolation

By definition, stress and force are connected by a geometrical area. In FEM, the conversion from stress to force is performed retrieving stress values at the Gauss points of mesh elements and converting them into nodal values using the pseudo-inverse of a shape function associated to each mesh type. Figure 3a shows a typical example about how this is implemented in El-Numodis, associated here to a a regular mesh of 8-nodes hexahedron elements for the sake of simplicity. As shown in Figure 3a, σ(S) is computed at the mesh nodes and the conversion from σ(S) to F (S) is done using a weighted area of normal vector n associated to each node using Equation 2. The weighted area depends on the element shape and localization. Here, internal surface nodes are characterised by A=1 whereas corners and external edge nodes have associated area of A/4 and A/2, respectively.

F (S) = Aσ(S)n (2)
The stress at the dislocation Gauss points σ gp is derived from the FEM solution σ originally computed at the mesh nodes. To interpolate and retrieve the solution at the dislocation Gauss point, it is necessary to (i) identify the element and the k nodes enclosing the respective dislocation Gauss point, (ii) transform the global coordinates (x,y,z) of the element nodes into a reference frame (ξ,η,ζ) (Figure 3b) and (iii) apply the shape function N k . In the square mesh case depicted in Figure 3a, N k and σ gp are provided by,

N k = 1 8 (1 + ξ gp ξ k )(1 + η gp η k )(1 + ζ gp ζ k ) (3)
σ gp = 8 k=1 N k σ k (4)
where k refers to nodes index. This method is commonly used in El-Numodis that benefits of various additional interpolation schemes (via Elmer) that could be adapted to the calculation of σ gp . See e.g., [START_REF] Zienkiewicz | The Finite Element Method: Its Basis and Fundamentals[END_REF] for more details on field conversion and interpolation.

Dislocation and surfaces

El-Numodis is designed to model dislocations interacting with surfaces (optionally) using the mirror dislocation concept [START_REF] Hirth | Theory of dislocations[END_REF][START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF] to assist the FEM in field corrections. Indeed, when a dislocation approaches a surface, the local stress and displacement fields can be obtained adding the contribution of an out-of-the-box mirror dislocation to the original dislocation fields. While the field calculation will be discussed later, the topological aspects of the mirroring process in El-Numodis can be described as follow. First, following the seminal work of Weygand et al. [START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF], a dislocation close to a sample surface by a cutoff distance r im c is automatically replicated on the other side of the surface using planar symmetry as shown in Figure 4. The resulting mirror dislocation is characterized by a line with symmetric orientation but an opposite Burgers vector direction. While the mirror dislocation stress and displacement fields are accounted within the simulation cell in order to reduce image contributions, the image dislocation does not produce any plastic shear. If the dislocation is about to contact with the surface, the Numodis collision detection algorithm is used to identify dislocation segments about to react with their mirrored counterpart, both emerging at the surface. In this case, the dislocation contact reaction leads to the annihilation of both dislocations. Finally, the dislocation annihilated portion at the surface is replaced by a ledge made of surface nodes (Figure 4). The surface nodes have the same mobility as the bulk ones but are constrained to move only on the sample surface (with the possibility to pass from one surface to an other). Also, ledges benefit from the same properties as dislocations i.e., they can superimpose or annihilate when several dislocations escape from the same surface but do not produce any stress or displacement field inside the sample.

Dislocation nucleation

Dislocation nucleation was originally introduced within DDD simulations by Fivel et al. who used a criterion on the macroscropic stress to model the nucleation of prismatic loops during nanoindentation [START_REF] Fivel | Three-dimensional modeling of indent-induced plastic zone at a mesoscale[END_REF]. More recently, Roy et al. investigated homogeneous dislocation nucleation in nanoparticles nucleating octagonal loops with characteristic size (imposed by the user) in high-stress regions further testing their expansion by a strict calculation of the Peach-Koehler force between two consecutive DDD time steps [START_REF] Roy | Annihilation of edge dislocation loops via climb during nanoindentation[END_REF]. This kind of approach was consecutively extended or adapted by several groups [START_REF] Gagel | Formation of extended prismatic dislocation structures under indentation[END_REF][START_REF] Po | Recent Progress in Discrete Dislocation Dynamics and Its Applications to Micro Plasticity[END_REF][START_REF] Roy | Annihilation of edge dislocation loops via climb during nanoindentation[END_REF]. At the atomic scale, SDN has shown to be a stochastic process that depends on both the local shear stress and temperature [START_REF] Zhu | Temperature and strain-rate dependence of surface dislocation nucleation[END_REF][START_REF] Li | Sample-size-dependent surface dislocation nucleation in nanoscale crystals[END_REF]. This kind of mechanism can be rationalized using the TST that connects the phenomenon frequency of occurrence to the local thermomechanical conditions through an activation energy, as it has been done e.g., by Zhu et al. for SDN in metal nano-objects [START_REF] Zhu | Temperature and strain-rate dependence of surface dislocation nucleation[END_REF]. Recently the approach of Zhu was extended to investigate the influence of the local atomic environment on the SDN process accounting for the sharpness of surface corners and edges [START_REF] Amodeo | Site dependence of surface dislocation nucleation in ceramic nanoparticles[END_REF][START_REF] Nisany | A Multiple Site Type Nucleation Model and Its Application to the Probabilistic Strength of Pd Nanowires[END_REF]. To provide a more realistic description of what is commonly done at the mesoscale, our approach benefits of the recent theoretical progresses made on the dislocation nucleation process. Indeed, the code uses a combination of dislocation harmonic TST and Kinetic Monte-Carlo (KMC) to identify the favorable homogeneous or heterogeneous dislocation nucleation sites using meshed activation data and a probabilistic approach. As suggested by the harmonic TST, the dislocation nucleation rate at finite temperature and site i is described by,

ν i = ν 0,i exp - ∆G i (σ i , T ) k B T (5)
where ν 0,i is the local attempt frequency, ∆G i is the Gibbs free energy of the dislocation nucleation process at site i, σ i is the local stress and k B T is the Boltzmann factor.

Thus, the probability for a dislocation to nucleate at site i during time δt in El-Numodis is computed as,

δN i = ν i .δt (6)
with

δt = 1 N i=1 ν i (7)
with N the total number of nucleation sites.

The frequency n KM C at which the KMC algorithm is called is defined by n KM C =δt DDD /δt per DDD step, where δt DDD is the DDD timestep. n KM C is adapted on-the-fly during the simulation based on possible ν i variations induced by local stress changes. Dislocation nucleation activation energy ∆G i and other characteristic parameters (critical nucleation radius and slip systems) can be set as non-local inputs or using tabulated data on a 2D or 3D grid to account for the site-dependence of the nucleation process. Corresponding dislocations are then generated into the simulation box as circular glissile loop (homogeneous nucleation) or truncated half or quarter loops (heterogenous nucleation) depending on the nucleation site location. Finally, El-Numodis is able to interpolate multiple stress or strain-dependent ∆G i databases to model the dislocation nucleation sensitivity to the mechanical history of the virtual sample.

While a comprehensive study of the SDN process in ceramic nanoparticles using atomistically-informed El-Numodis will be the main focus of a forthcoming study, a simplified application of the dislocation nucleation process is presented in Figure 5.

In this example, we use a cubic-shaped copper single crystal of 500×500×500 nm 3 size under constant load. SDN in the 1 2 110 {111} slip systems at three hypothetical nucleation sites is considered i.e., two top corners c 1 and c 2 and the middle of a single lateral surface labelled s. As in the work of Zhu et al. [START_REF] Zhu | Temperature and strain-rate dependence of surface dislocation nucleation[END_REF], we use the approximation of an homogeneous surface disordering temperature ∆S i = ∆H i /T m,i , where ∆S i is the activation entropy at site i, ∆H i is the 0 K nucleation activation enthalpy and T m,i is the local surface disordering temperature (ν 0,i =3.14 10 11 /s and T m,i =700 K whatever i). Corner and mid-surface activation enthalpies are set to ∆H c 1 =∆H c 2 =0.2 eV and ∆H s =0.5 eV respectively, assuming SDN to be more efficient from corners than from mid-surface. A high-enough constant load is considered (σ zz =1 GPa) to guarantee dislocation nucleation and glide. For a sake of simplicity, the per-site activation nucleation energy is temporarily increased after each nucleation event to mimic the effect of internal stress relaxation that prevents overly correlated nucleation events. Simulation performed at T =250 and 650 K are described Figure 5. El-Numodis promotes SDN originating only from c 1 and c 2 at T =250 K (Figure 5a) while nucleation events incoming from the three sites are observed at larger temperature (Figure 5b). This result is mainly explained by the site-dependence of the nucleation probability that is influenced by the temperature range as expected by the harmonic TST. Here, the probability to nucleate from a corner P c =ν c / i ν i ∼0.5 at T =250 K while its mid-surface counterpart is close to zero (P s =6. 47 10 -5 ) what justifies the lack of nucleation event from site s at low temperature. On the other hand, P s increases up to ∼0.25 when the temperature is increased up to 650 K as shown Figure 5c and nucleation from site s becomes more favorable.

Loading and feedback loops

El-Numodis handles displacement-and force-controlled BCs (Dirichlet and Neumann BCs, respectively) that can be set at each surface of the virtual sample allowing for constant strain rate or creep simulations. During constant strain rate simulation, a feedback control acting on stress is performed as described in Equation 8 using Neumann BCs. On the other hand, the sole SPM method is used during creep simulations.

σ(S, t) = 1 Σ ( ε.t -ε p (t)) (8)
where σ(S, t) is the applied stress tensor at the surface, t is the elapsed time, ε is the total strain rate tensor and ε p is the plastic strain tensor computed using the area swept by all moving dislocations. Σ is the corresponding compliance tensor.

Validation and applications

Dislocation stress field

El-Numodis benefits of both the singular and non-singular formulations for dislocation stress field (equations are provided as Supplementary Information). Figure 6 shows a comparison between El-Numodis DDD and the analytical formulation of the stress field using the Cai's non-singular theory [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF].From a technical point of view, a straight (edge or screw) dislocation is modeled assuming an infinite continuum using a simulation cell of size 1×1×1 µm 3 . The dislocation line is oriented along z=[001] and the stress field is illustrated in the (x,y) plan at z=0.5 µm. Copper lattice parameter (a 0 =3.61 Å) and isotropic elasticity (λ=77.3 GPa, µ=42.0 GPa, ν=0.324) are used for the example. Overall, results show a good agreement between analytical and El-Numodis solutions. One could notice that changing the BCs from periodic to fixed-BCs as well as short variations of cell dimensions along x and y do not significantly impact the results for the investigated size range.

Edge dislocation and free surface

Here we test El-Numodis reliability in the context of the stress-free BC problem by investigating how the stress field of an edge dislocation is modified in the vicinity of a free surface, as described in Figure 7a. In their seminal work, Hirth and Lothe proposed an analytical solution to this problem that relies on the concept of image dislocation [START_REF] Hirth | Theory of dislocations[END_REF]. Indeed, the authors demonstrate that most of the stress components induced by an infinite edge dislocation can be cancelled at a surface by adding a socalled image dislocation located on the other side of the surface. The image dislocation is characterized by i) the same line direction and infinite length, ii) an opposite Burgers vectors -b and iii) the same dislocation-to-surface distance l than the original dislocation. In addition, they use an Airy stress σAiry superimposed to the self σ and image dislocation σim stress-fields to ensure the complete vanishing of all the stress components at the free surface (i.e., for x=l), without modifying the long-range stress distribution within the sample,

σij + σim ij + σAiry ij = 0, for x = l (9)
with σim ij the stress components of the image dislocation.

Within the Hirth and Lothe 2D formulation, the Airy stress components that verify the stress-free conditions of the x-oriented surface are given by,

σ Airy xx = - 2µblxy π(1 -ν)r 6 [3(l -x) 2 -y 2 ] ( 10
)
σ Airy xy = - µbl π(1 -ν)r 6 [(l -x) 4 + 2x(l -x) 3 + 6xy 2 (l -x) -y 4 ] (11)
where r = (l 2 + y 2 ) 1/2 The Hirth and Lothe model presented in Equation 9 relies on a stress summation very close the one used in the SPM. Thus, one way to test El-Numodis implementation is to verify that the FEM stress correction computed by the code correctly reproduces the Hirth and Lothe theoretical predictions of Airy and image stresses i.e., σ=σ im +σ Airy . To test this hypothesis, we design 1×2×1 µm 3 simulation cells including a finite-length edge dislocation (line along z, Burgers vector along x) located at a distances l=0.1 µm from a surface located at x=0. Simulation volumes are meshed using 8-nodes hexahedron elements. Various geometries are tested including elements distributions from 30×30×30 up to 120×120×120. Figure 7a shows a 30×60×30 volume with mesh refinement near the dislocation performed using Gmsh (Bump=-4.9 and Progression=1.3) [START_REF] Geuzaine | Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. As in the previous test-case, copper lattice and elastic properties are used. In the simulation, traction-free BC is applied to the x free-surface while zero fixed displacement is used on the opposite side (other surfaces are not considered). The total stress σ and the dislocation self-stress σ obtained are presented in Figure 7b as plotted along a vertical line passing along the free-surface. Results confirm El-Numodis ability to reproduce theoretical σ, as already shown in Figure 6. In addition, the plot of σ allows for a direct comparison between the FEM correction σ and the Hirth and Lothe model σim +σ Airy . The total stress components σ xx and σ xy show a significant decrease at the surface (originally equal to the self-stress) what confirms the correct implementation of the superposition algorithm within El-Numodis. This result applies almost everywhere except close to y=0, where spurious stresses are observed. These discrepancies are attributed to the mesh refinement in this region (where large stress gradients are shown) as well as to the intrinsic difference between the 2D Hirth and Lothe model and the 3D finite-size simulation. One can notice that similar stress singularities at free-surface were already observed in simulations using the SPM [START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF] or the DCM approach [START_REF] Jamond | Consistent formulation for the Discrete-Continuous Model: Improving complex dislocation dynamics simulations[END_REF][START_REF] Huang | An extended 3D discrete-continuous model and its application on single-and bi-crystal micropillars[END_REF].

Still for this test-case, Figure 8 illustrates the influence of the mirror dislocation method on El-Numodis stress field correction. Results show a qualitative agreement between simulated σ (or σ-σ im , depending if the mirror image method is turned off or on) and the σAiry of the Hirth and Lothe model. Overall, using the mirror image method with FEM to compute the image stress improves the results allowing for less refined meshes. However, such an improvement has limitations as emphasized by the y = 0 region where data for 60 and 120 elements simulations saturate. Thus, the residual stress in the y = 0 region also observed Figure 7b is neither significantly sensitive to simulation cell size variations or nor to the mesh refinement.

Square dislocation loop and free surface

In this section, we use El-Numodis to solve the stress-free BC problem in the case of a dislocation square loop located in the vicinity of a free surface. Curved dislocation stress field can be computed using linear elasticity and surface integration built out of the dislocation curvature [START_REF] Fivel | Implementing image stresses in a 3D dislocation simulation[END_REF][START_REF] Mura | Continuous distribution of moving dislocations[END_REF][START_REF] Gosling | A line-integral representation for the stresses due to an arbitrary dislocation in an isotropic half-space[END_REF],

σ pq (x ) = -b s C srkl C pqmj ∂ ∂x j G mk,j (x, x )dS r ( 12
)
where C is the elastic stiffness and G is the Green tensor associated to a particular material.

In an infinite medium, Equation 12 can be reduced into a simple integral computed along the dislocation line using the Stokes theorem. Gosling and Willis expanded this approach to a finite-size domain using Equation 13where S ∞ and Ŝ are kernels associated to infinite and finite-size media respectively [START_REF] Gosling | A line-integral representation for the stresses due to an arbitrary dislocation in an isotropic half-space[END_REF]. The integral of Ŝ directly leads to the image stress σ; It will be referred as the Gosling-Willis solution in the following (more details on this approach can be found e.g., in Ref. [START_REF] Fivel | Implementing image stresses in a 3D dislocation simulation[END_REF]).

σ pq (x ) = - C b s [S ∞ pqrs (x -x ) + Ŝpqrs (x, x )]dx (13)
The El-Numodis simulation setup used to model the square loop and surface interactions is presented Figure 9a. A copper 100 -oriented simulation cell with dimensions [ 101](111) square dislocation loop axis-aligned with z and with edge lengths of about 0.5 µm is introduced at 0.37 µm of the bottom-z surface at which stress-free BCs are applied. As in the previous case, zero displacement fixed BCs are used for the opposite z surface while remaining surfaces are not considered as boundaries in the simulation. The mirroring dislocation method is off here.

Results are presented in Figure 9b and 10. On one hand, Figure 9b shows the variations of the FEM stress correction σ as plotted along the z-axis starting from the free surface up to 2 µm and passing trough the center of the dislocation loop. Here the results illustrate that El-Numodis is particularly suited to reproduce the Gosling-Willis theoretical solution. On the other hand, σ, σ and σ+σ stress maps plotted at the bottom-z surface are shown in Figure 10. El-Numodis alllows for the decrease of the various stress components by more than a factor 5 with only few MPa leftover at the surface. This test further confirms the ability of El-Numodis to solve stress-free BVPs.

Example of application: thin film tensile test

For this final application, El-Numodis is confronted to classical DDD (i.e., standalone Numodis) performing tensile tests on a model thin film. A 100 -oriented Cu thin film of 250 nm thickness is generated and meshed using eight nodes hexahedrons discretized with 20 points along [010] and [001] directions while 10 points are used along [100]. The initial dislocation microstructure is identical in both simulations. It is made of 25 Frank-Read sources of 0.2 µm length randomly distributed on the various slip systems of the FCC crystal structure. The defects distribution is biased using a cutoff to avoid overlaping and boundary crossing. Constant strain-rate simulation ( ˙ =10 -6 /ns) is performed pulling from one of the (001) lateral surface while the opposite one is kept fixed. Other surfaces are set stress-free in El-Numodis. The mirror-image dislocation method is used in El-Numodis simulation only with a cutoff distance of 60 nm while free-BCs are used in the pure DDD simulation. Equation 8 feedback loop is used to correct the applied stress at the boundary in both simulations. Computed stress-strain curves are shown in Figure 11a. Both SPM and DDD simulations are characterized by an initial elastic load up to the activation of the first Frank-Read sources. The pure DDD exhibits a higher yield stress when compared to El-Numodis (196 and 175 MPa, respectively) and, overall, a harder mechanical response Here surfaces in the El-Numodis simulation behave as sinks that help the opening of the Frank-Read sources towards the surfaces in a similar way that in the aforementioned TEM lamella case [START_REF] Aaron A Kohnert | On the use of transmission electron microscopy to quantify dislocation densities in bulk metals[END_REF]. Moreover, accounting for surfaces significantly influences the dislocation dynamics and elementary reaction processes between dislocations. For example, some junction reactions initially observed in the pure DDD simulation have shown to be anticipated by the surface-induced pump-out process modeled using El-Numodis and the SPM method. This process influences the whole plastic regime which is shown to be softer when accounting for the physics of surfaces. Additional simulations show that increasing the film thickness tends to reduce the gap between pure DDD and the SPM mechanical responses while increasing the cutoff used in the mirroring dislocation method tends to decrease the yield stress.

Conclusion

Here we present a tool based on the SPM method called El-Numodis to model dislocation dynamics in finite-size environments. The approach couples the nodal DDD code Numodis and the Elmer FEM code used here as an elastic solver of BVPs where the BCs are corrected by the effect originating from the presence of dislocations. El-Numodis refers on three external drivers that ensure the various operations and exchange between the DDD and FEM parent codes. It benefits of specific modern developments including the non-singular dislocation theory of Cai [START_REF] Cai | A non-singular continuum theory of dislocations[END_REF], the mirror image method [START_REF] Weygand | Aspects of boundary-value problem solutions with three-dimensional dislocation dynamics[END_REF] as well as a Monte-Carlo based dislocation nucleation algorithm that allow for more physics-based dislocation simulations at the micro-and nano-scales. From a technical viewpoint, these developments make El-Numodis particularly versatile (and not more complex) than its original DDD parent code Numodis [START_REF] Drouet | Dislocation dynamics simulations of interactions between gliding dislocations and radiation induced prismatic loops in zirconium[END_REF][START_REF] Shi | Interaction of < 100 > dislocation loops with dislocations studied by dislocation dynamics in α-iron[END_REF]. In this study, El-Numodis was widely benchmarked including several test-cases designed to investigate interactions between dislocation and surfaces. Among others, straight and square-shaped dislocations were tested in the vicinity of free-surfaces where El-Numodis has shown to be particularly suited to relax surface stress fields. A dislocation nucleation algorithm using a Monte-Carlo approach and the TST was also introduced. It will be soon extended to the study of dislocation nucleation in nanoparticles using energy barrier databases [START_REF] Amodeo | Site dependence of surface dislocation nucleation in ceramic nanoparticles[END_REF]. Finally, last tensile tests applications performed on model thin films completes the picture of El-Numodis potential. Besides ongoing improvements as on force calculation [START_REF] Queyreau | Analytical integration of the forces induced by dislocations on a surface element[END_REF][START_REF] Queyreau | Analytical integration of the tractions induced by non-singular dislocations on an arbitrary shaped triangular quadratic element[END_REF] or Gauss integration [START_REF] El-Achkar | Aspects on numerical integration of dislocation surface traction fields for discrete dislocation dynamics fem coupling: the case of emerging dislocations[END_REF], we believe that El-Numodis is now ready for various kinds of applications in the field of small-scale mechanics.

Figure 1 :

 1 Figure 1: The superposition method. The stress field σ inside a finite-size volume containing dislocations (right-hand side) is obtained adding the dislocation self-stress σ as computed by the DDD to the elastic stress σ as inferred from a FEM boundary corrected problem (left-hand side).

Figure 2 :

 2 Figure 2: El-Numodis workflow. The FEM code Elmer masters the coupling using the DDD Numodis code as a library. Both contributions are illustrated in green and blue, respectively. The three main drivers of El-Numodis i.e., El-Numodis export Bnodes, El-Numodis import and El-Numodis export stress are illustrated using red blocks. (i), (ii), (iii) and (iv) refer to the main operations as described in the text.

Figure 3 :

 3 Figure 3: El-Numodis interpolation scheme using shape functions and fields conversion. (a) Representation of a surface mesh where stresses are converted into forces at the element nodes. Nodes and surface colors refer to respective types (green=corner, blue=side edge and red=internal). (b) The global coordinates (x,y,z) of each node is transformed into reference coordinate (ξ,η,ζ) with the new axes located at the center of the hexahedron element.

Figure 4 :

 4 Figure 4: Dislocation interacting with surfaces in El-Numodis. (a) a dislocation close to a surface by a distance r im c is mirrored using planar symmetry, (b) bulk and mirror dislocations react at the surface creating a surface ledge. Dislocation and ledge nodes are colored according to their Degree Of Freedom number (DOF). NA=Not Accounted.

Figure 5 :

 5 Figure 5: Dislocation nucleation at the surfaces of a cuboïdal Cu sample under constant load (σ zz =1 GPa) at (a) T =250 K and (b) T =650 K. (c) Nucleation probability P i =ν i / i

Figure 6 :

 6 Figure 6: Comparison between El-Numodis and analytical formulations of edge and screw dislocation stress fields in the framework of the non-singular theory of Cai [50].

Figure 7 :

 7 Figure 7: Edge dislocation and free surface: stress-field prediction comparison between El-Numodis simulation (Sim) and Hirth and Lothe (H&L) theoretical model. (a) Simulation setup, (b) Dislocation self-stress σ and total stress σ variations along y for a 30×60×30 elements mesh with special refine in the vicinity of the dislocation (Gmsh Bump=-4.9 and Progression=1.3). Data are plotted along a virtual line crossing the free surface.

Figure 8 :

 8 Figure 8: El-Numodis σxy (or σxy -σ im xy when accounting the mirroring dislocation method) correction as function of mesh refinement. Simulations are performed for mesh discretization of 30 (dashed curves), 60 (dotted curves) and 120 (full curves) 8nodes hexahedron elements in the three directions of space, using the Weygand's mirror dislocation method (W) or not. Ref curves rely on the aforementioned 30×60×30 with particular mesh refinement near the dislocation using Gmsh (Bump=-4.9 and Progression=1.3). Results are compared to the σAiry xy of the Hirth & Lothe model (H&L, red curve). Data are plotted along a vertical line passing by the middle of the x-surface as shown in Figure 7a.

Figure 9 :

 9 Figure 9: Dislocation square loop in the vicinity of a free surface: image stress field comparison between El-Numodis simulation (Sim) and the Gosling-Willis (GW) theoretical model (a) Simulation setup, (b) Image stress components σ ij as measured along a vertical line passing by the middle of the dislocation square loop. El-Numodis results (dots) are compared to the Gosling-Willis model (curves).

Figure 10 :

 10 Figure 10: Stress components of σ, σ and total stress σ+ σ computed using El-Numodis at the bottom-z surface close to a dislocation square loop.

Figure 11 :

 11 Figure 11: Thin film tensile test: standalone Numodis (DDD) vs. El-Numodis (SPM). (a) stress-strain curves, (b) Comparative dislocation microstructures at the yield point (black=Numodis, green=El-Numodis). Red arrows highligth early emerging dislocations in the El-Numodis simulation. Black arrows refer to the tensile direction.

Acknowledgements

This work was founded by INSA-Lyon and the Agence National de Recherche, grant no. ANR-20-CE09-0015 (ANR SASHA).

Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files). El-Numodis is available on-demand at laurent.dupuy@cea.fr and jonathan.amodeo@cnrs.fr.


Arthur Bongrand

Samuel Duval

Jérôme Donnard

Julie Champion

Latest development of

Moreover, the spatial resolutions of the two methods are often detector dependent, which makes it difficult to integrate and interpret these two measurements together.

Can we perform a fast and accurate quality control of radiopharmaceuticals with a single system?

References

R.O.I. 1D projection

In nuclear decay physics, a radioactive equilibrium exists when a radioactive nuclide decays at the same rate as it is produced.

Bateman equation (general explicit formula)

→ The concentration of daughter nuclei in the radioactive equilibrium can be modeled at any time using the Bateman equation provided that the characteristics of the decay chain and if the initial concentrations are known. → For our application, as these initial concentrations are unknown, we have developed a method to calculate these concentrations from the evolution of the measured count rate.

Count mode

Mean Time mode

Mouse kidney labelled with 225 Ac-PSMA

→ Usually the quality controls are made from a 1D projection where we can distinguish the free 225 Ac which remains on the deposit spot from the labelled 225 Ac which migrates. To limit interpretation errors, the analysis must be performed at secular equilibrium (~7h).

→ The STSAA reduces the time needed for a quality control after elution to only one hour allowing both isotope localization and identification. Thus, the purity is deduced from the ratio of free to labeled 225 Ac.

Digital autoradiography combined with a temporal spectroscopy method allows the quick and accurate identification of different radioisotopes and to observe their localization.

The increased robustness of this method comes from its low dependence on the operator's choices compared to 1D method.

From the results, we can observe the contribution of a radioisotope called "free 213 Bi". Clearly visible on the migration front, free 213 Bi corresponds to 213 Bi which is not radiolabeled and whose parent, 217 At, has not decayed locally. The most probable hypothesis is that the free 213 Bi comes from the mother solution. This free 213 Bi has also been observed in biological samples. Studies are underway to understand its biological impact.

Results

Discussion

  𝑖 𝑁 𝑖 𝑡 + 𝜆 𝑖-1 𝑁 𝑖-1 𝑡 𝑑𝑁 𝑘 𝑡 𝑑𝑡 = 𝜆 𝑘-1 𝑁 𝑘-1 𝑡 International Symposium on Trends in Radiopharmaceuticals (ISTR-2023), CN-310 Vienna, Austria; 17-21 April 2023 Latest development of α emitter imaging and quantification on a large Field Of View A. BONGRAND 1 , S. DUVAL 2 , J. DONNARD 2 , J. CHAMPION 1 1: Laboratoire Subatech, IMT Atlantique, Nantes Université, CNRS/IN2P3, France 2: S.A.S. Ai4R, France arthur.bongrand@subatech.in2p3.fr Currently, the characterization of the nature and spatial distribution of radionuclides in a sample is time consuming and fastidious. Indeed, it necessarily requires the use of two distinct analytical techniques and detectors of different nature.

✓✓

  Promising in the acquisition of nuclear data for the production of radionuclides for nuclear medicine applications like quality control of radiopharmaceuticals ! ✓ A precise tool to understand the temporal behavior of a chain of labeled isotopes. Precise image & quantification from mBq to kBq. More experiments ! → Perform simultaneous comparisons with standard quality control methods. → Perform multiple quality controls of the same batch simultaneously using the BeaQuant's large Field Of View. → Explore the method for other radiopharmaceuticals [2]. up-to-date emission imaging technique ! ✓ Allows to study and locate in 2 dimensions the radioactivity emanating from a surface sample at scales until 10 micrometer. ✓ Well suited for pre-clinical TAT as it allows understanding and quantifying biodistributions in tissues at high spatial resolution and low activity. ✓ Digital autoradiography → Ability to analyze and reanalyze endlessly the digital data acquired! ✓ Variety of applications: Radio-TLC, Geological Studies, Radio Toxicology…

Acknowledgments

We thank the GIP Arronax for providing of their facilities and a 225 Ac-PSMA source, especially Dr. Cyrille Alliot for his excellent expertise in the field of radiopharmaceutical quality control.