Rh(III) and Ir(III)

Subscriber access provided by UNIV OF NEWCASTLE

Article

DFT Simulation of Structural and Optical Properties of 9Aminoacridine Half-Sandwich Ru(II), Rh(III) and Ir(III) Antitumoural Complexes and their Interaction with DNA José Pedro Cerón-Carrasco, José Ruiz, Consuelo Vicente, Concepcion de Haro, Delia Bautista, Jose Zuniga, and Alberto Requena J. Chem. Theory Comput., Just Accepted Manuscript • Publication Date (Web): 22 Jun 2017 Downloaded from http://pubs.acs.org on June 23, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

DFT Simulation of Structural and Optical Properties of 9-Aminoacridine Half-Sandwich Ru(II), Rh(III) and Ir(III) Antitumoural Complexes and their Interaction with DNA José Pedro Cerón-Carrasco,∗,† José Ruiz,‡ Consuelo Vicente,‡ Concepción de Haro,‡ Delia Bautista,¶ José Zún˜iga,§ and Alberto Requena§ †Bioinformatics and High Performance Computing Research Group (BIO-HPC) Universidad Cat´ olica San Antonio de Murcia (UCAM) Campus de los Jerónimos, 30107, Murcia, Spain ‡Departamento de Qu´ımica Inorg´ anica, Universidad de Murcia, 30100 Murcia, Spain ¶SAI, Universidad de Murcia, 30100 Murcia, Spain §Departamento de Qu´ımica F´ısica, Universidad de Murcia, 30100 Murcia, Spain E-mail: [email protected]

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract In this work, we use DFT-based methods to simulate the chemical structures, optical properties and interaction with DNA of a recently synthesized chelated Cˆ N 9-aminoacridine arene Ru (II) anticancer agent and two new closely related Rh(III) and Ir(III) complexes using DFT-based methods. Four chemical models and a number of theoretical approaches, which representatively include the PBE0, B97D, ωB97X, ωB97X-D, M06, and M06-L density functionals and the LANL2DZ, def2-SVP, def2TZVP basis sets, are tested. The best overall accuracy/cost performance for the optimization process is reached at the ωB97X-D/def2-SVP and M06/def2-SVP levels of theory. Inclusion of explicit solvent molecules (CHCl3 ) further refines the geometry, while taking into account the crystal network gives no significant improvements of the computed bond distances and angles. The analysis of the excited states reveals that the M06 level matches better the experimental absorption spectra, compared to ωB97X-D. The use of the M06/def2-SVP approach is therefore a well-balanced method to study theoretically the bioactivity of this type of antitumoral complexes, so we couple this TD-DFT approach to molecular dynamics simulations in order to assess their reactivity with DNA. The reported results demonstrate that these drugs could be used to inject electrons into DNA, which might broaden their applications in photoactivated chemotherapy and as new materials for DNA-based electrochemical nanodevices.

Keywords: Organometallics, DNA intercalators, anticancer compounds, density functional theory, benchmark.

2


Page 2 of 39

Page 3 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

The increasing development of the density functional theory (DFT) methods means it is now possible to support the experimental study of the geometry, optical signatures and chemical reactivities of newly synthesized complexes by conducting accurate computational simulations with the aim of anticipating the anticancer activity of promising derivatives. 1–6 In this respect, we previously used a wide panel of theoretical approaches to model the cytotoxicity of platinum-based drugs. 7–9 It is also widely accepted that the ultimate step in the antitumoral properties of metallodrugs is the DNA. 10–13 However, despite all the accumulated data on the organometallic complexes’ reactivity with DNA, the experimental techniques available are not able to resolve the reaction mechanisms. For this purpose one needs a robust methodology to accurately mimic the interaction of the metallodrug with the biological molecules (i.e., DNA or protein), which presumably requires the combination of different levels of theory, including molecular dynamics (MD), time-dependent DFT and hybrid quantum mechanical-molecular mechanics (QM/MM) and quantum-mechanical-semi empirical (QM/SE) calculations. 7,14,15 The accuracy of a given theoretical approach is closely related to the chemical system under study, so the computational protocol applied must be carefully selected. For organometallic derivatives, the weak nature of the metal–carbon bond is a major handicap for the calculation of the electronic properties. 16 Since no DFT method can be considered reliable for all metal derivatives, 17,18 the accuracy of the theoretical approach has to be checked in each specific case in order to make valuable predictions. Ruthenium derivatives present many possible applications in a wide panel of fields, including biochemistry, 19 electrochemical materials, 20 charge-transfer sensitisers 21 and catalysis. 22 In the framework of novel organometallic compounds with biological applications, we have recently experimentally designed a novel Ru(II) intercalator which exerts a potent growth inhibitory effect on tumor cells (Complex-I in Figure 1). 23 In addition, two new closely related Rh(III) and Ir(III) derivatives (Complex-II and Complex-III in Figure 1, respectively) 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

C24

C21 Ru C1

C21

C24

Cl

Rh

N9 H

C1

C9

Cl N9 H

Page 4 of 39

C21

C24 Ir

C1

N

H Complex-I

H Complex-II

N9 H C9

C9

N

Cl

N H Complex-III

Figure 1: Chemical structures of the ruthenium compounds discussed in the text. have been ad-hoc synthesized to cover other isoelectronic metallic centers, and their X-ray crystal structures have been determined. In this work, we carry out a critical assessment of the efficiency of several functionals, basis sets, and chemical models for describing the geometrical parameters and absorption spectra of the 9-aminoacridine arene Ru(II), Rh(III) and Ir(II) complexes. As we are mainly focusing on the potential anticancer activity of these novel intercalators, it is also interesting to determine whether the novel drug can be photoactivated once the malignant DNA has been reached. 24–28 Accordingly, we conduct additional MD simulations in order to determine DNA binding modes of the selected organometallic compounds. Finally, TD-DFT 29 calculations are carried out in the DNA-embedded model systems to assess the capacity of the theory to describe the excitation processes and their application as photoactivatable drugs.

2

Chemical models and computational approaches

We have focused on Complex-I in order to evaluate the impact of the DFT approaches in the geometry optimization. This complex crystallizes with two solvent molecules (CHCl3 ) with their hydrogen atoms orientated towards the Cl− counter-ion, and shows dimeric units, likely due to a π-stacking interaction between two adjacent aminoacridine groups. 23 This suggests 4


Page 5 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


!"#$%&!

!"#$%&"

!"#$%&#

!"#$%&$

Figure 2: Selected model systems to determine the crystal environment effects, shown as balls and sticks representation. Model D, which is treated with the hybrid QM/MM and QM/SE approaches is given as balls and sticks (high layer) and wireframes (low layer). See the text for further details. that the complex is highly interconnected with the crystal surrounding environment, making it an excellent candidate for testing the importance of the chemical models on the balance between experimental data and theoretical predictions. Accordingly, we have built up four chemical models of increasing complexity for ComplexI from its X-Ray structure, which are shown in Figure 2. Model A corresponds to the isolated Complex-I structure (dry monomer); model B adds the two CHCl3 molecules interacting with the Cl− counter-ion (solvated monomer); model C represents a dimeric unit with the four solvent molecules shell (solvated dimer); and model D includes the solvated dimer in contact with 12 solvated monomers (12 metallic centres) accounting for the crystal network. The four model systems have been designed using Mercury CSD. 30 We have used the smallest model A to test the level of theory in the structure optimization. The geometry of the isolated Complex-I has been fully optimized, in particular, by combining the PBE0, 31 B97-D, 32 ωB97X, 33 ωD97X-D, 34 M06, 35 M06-L 36 density functionals with the LANL2DZ, 37 def2-SVP and def2-TZVP basis sets. 38 The choice of these functionals and basis sets is guided by previous benchmarks performed on Ru(IV) complexes. 39–43 The functionals selected cover pure, hybrid and range-separated hybrid schemes, that also allow us to assess the impact of the dispersion correction on the performance of the density func5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

tionals, 44 whereas the basis sets proposed allow us to explore size effects by comparing the results obtained through double-ζ valence (DZ) [also called split valence (SP)] and triple-ζ (TZ) valence basis sets. The core electrons of heavy atoms are replaced by effective core potentials (ECPs) to reduce the number of basis functions, while simultaneously accounting for scalar relativistic effects. The optimized structure is confirmed as a true minimum in the potential energy surface through the vibrational analysis of the stationary point (no imaginary frequencies). The rest of the models proposed can also be used to assess the influence of the method in the optimization procedure, but at a higher computational cost (larger systems). Accordingly, models B and C are optimized only with the benchmarked method, and model D is partially optimized using a hybrid (ONIOM) approach 45 due to its much larger size (770 atoms). This constrained ONIOM optimization has been previously shown by Kami´ nski textitet al. 46 to be a suitable procedure to explore the structure of a metallic derivate embedded in a crystal. On applying it to model D, we split the system into two regions or layers, a central region containing the solvated ruthenium dimer (high layer), which is fully relaxed using the benchmarked method, and an external region formed by the rest of the crystal environment (low layer), which is frozen in space. The low layer is described using molecular mechanics (MM) with the Universal Force Field (UFF), 47 and also the semiempirical (SE) PM6 method, 48 through their corresponding QM/MM and QM/SE strategies, respectively, with the former QM/MM electronic embedding scheme incorporating the partial MM charges in an additional term of the QM Hamiltonian. 49 To ascertain the accuracy of that hybrid scheme, model D is also fully optimized at the QM/MM level. In order to understand the chemical structure and photochemical properties of the three complexes, we have used TD-DFT to simulate their absorption UV-vis spectra. The TDDFT calculations were conducted using the density functionals that produce the most accurate structures, and accounting for the solvent effects with the Polarisable Continuum Model (PCM) in its linear-response non-equilibrium regime. 50 The nature of the excited

6


Page 6 of 39

Page 7 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


!"#$%&'('#)*"+'"-+&*0'($"#+,)"-)".+/*-$

!"#$%&'('#)*"+,)"-)".+/*-$

Figure 3: Representation of the chemical systems used in this study. Middle panel shows the initial DNA–drug adducts, where DNA is represented by ribbons and the intercalated Complex-I is drawn with blue and red balls-and-sticks, representing pure intercalation and combined intercalation and covalent binding modes, respectively. Left and Right panels zoom into the intercalation site. states involved in the main optical signatures were characterized through the charge transfer (CT) parameters defined by Le Bahers and co-workers. 51,52 In this analysis, the electronic transitions are represented by the differences between the total densities of the excited and ground states, which allow computation of the amount of charge transferred (q CT in |e|) and the distance associated to the electronic rearrangement (dCT in ˚ A). 51,52 All the geometry optimisation and absorption spectra calculations were carried out using Gaussian 09. 53 Since the ultimate goal of the computational simulations is to predict the anticancer activity of the metallodrug prior to its synthesis and in vivo/in vitro tests, it is necessary to design an accurate model system to mimic both the stability and the photochemical properties of DNA–drug adducts under biological conditions. In the present case, an initial DNA system was formed starting from the experimental X-ray crystal structure with the d(CGCGAATTCGCG)2 sequence as deposited in the Protein Data Bank, code 1G3X. 54 We chose this dodecamer because it advantageously includes an acridine-based drug (actinomycin D) intercalated in the middle of the base pairs, so the initial double helix architecture is ready 7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to host the three-member rings of the acridine moiety when conducting the MD simulations. In addition, the central intercalation avoids end-effects arising from the border base pairs. 9 Middle panel of Figure 3 shows that there are two possible DNA–drug binding modes: a pure non-covalent intercalation, in which the acridine moiety docks between two consecutive DNA basis pairs, and a combined intercalation/covalent binding. In the former, the intercalator present in the original PDB file, actinomycin D, is replaced by Complex-I by matching both acridine ligands. It is well established that inside the cell a water molecule is incorporated into the complex by displacing the chlorine ligand. 55–58 As shown in the left panel zoom, the water molecule of the resulting activated aqua-complex is retained in the intercalation binding mode. As for the combined intercalation/covalent binding mode, it arises from the common reactivity of metallodrugs, according to which they may react with the N7 center of the guanine (G) in the guanine-cytosine (GC) base pair. 59–61 The right panel of Figure 3 shows in this case that the during DNA attack, the position of the water molecule is exchanged by the N7 atom. Accordingly, the combined intercalation/covalent binding mode is built up by replacing the water molecule linked to the metal by the N7 of the vicinal guanine base. As illustrated in the middle panel of Figure 3, the reaction with the N7 position (drug in red) requires a slight twist of the acridine ring compared to the aqua-complex (drug in blue). Both binding modes were assessed by locating the two possible adducts in an orthorhombic box, with a buffer distance of 10 ˚ A in all directions, and by subsequently filling the box with water molecules SPC model and sodium cations to keep the system electronically balanced. Additional sodium and chloride ions were incorporated into the system to simulate the physiological salt concentration of 0.15 M NaCl. As far as the MD simulations is concerned, an initial full minimization was conducted for 2000 steps using the steepest descent method, with a convergence threshold of 1.0 kcal/mol/˚ A. The solvated systems were next relaxed throught several stages that include a solute-restrained minimization, free-restrain minimization, NVT simulation of 24 ps at T=10 K, and NPT simulations at T=10 K and P=1 atm. For the production phase, the

8


Page 8 of 39

Page 9 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


temperature was set to 300 K with the Nosé-Hoover algorithm, with a relaxation time of 1.0 ps. 62,63 Pressure was controlled to 1 bar with the Martyna–Tobias–Klein barostat using isotropic coupling and a relaxation time of 2.0 ps. 64 The RESPA integrator was used to integrate the motion equations with a 2.0 fs time step for bonded and near interactions, and a 6.0 fs time step for far interactions. 65 A cut-off of 9 ˚ A was applied to non-bonded interactions. Van der Waals interactions were evaluated using a cut-off radius of 9 ˚ A, and the electrostatic part was computed using the Particle Mesh Ewald (PME) 66 method with a tolerance value of 10−9 . MD simulations were run for 20 ns using the recently developed OPLS3 force field as implemented in Desmond code, 67,68 which has been successfully used for simulating DNA intercalators. 69 Since no parameters are available for these new transition metal complexes, bond and angle constraints were imposed to the metal coordination sphere during the MD production. In order to relax the restricted metal environment, the chemical region of interest, the sandwiched Complex-I and the border basis pairs, were eventually extracted from the last MD snapshot to be further refined with the benchmarked QM method, and to determine the photochemical properties of damaged DNA by using the benchmarked TDDFT protocol. Finally, we applied the benchmarked TD-DFT protocol to the DNA–drug model to investigate the possible application of the present triad of metallodrugs in the photoactivated chemotherapy (PACT) framework, which offers an interesting alternative to partially control the drug activation process.

3

Geometry optimization

Let us consider first the geometry optimisation of Complex-I, as carried out comparatively using all the levels of theory described above. We select for this purpose the bond lengths and angles employed to characterize the experimental X-Ray structure of the complex, 23 namely, the Ru–Cl, Ru–C1 , Ru–C21 , Ru–C24 , Ru–N9 and N9 –C9 bonds and the C1 –Ru–N9 ,

9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 39

Table 1: Unsigned Error for Selected Bond Lengths and Mean Unsigned Error (MUE) for Optimized Model A (in ˚ A)a bond Ru–Cl

basis set

PBE0

LANL2DZ 0.009 def2-SVP 0.057 def2-TZVP 0.063 Ru–C1 LANL2DZ 0.040 def2-SVP 0.027 def2-TZVP 0.025 Ru–C21 LANL2DZ 0.080 def2-SVP 0.002 def2-TZVP 0.006 Ru–C24 LANL2DZ 0.119 def2-SVP 0.010 def2-TZVP 0.006 Ru–N9 LANL2DZ 0.035 def2-SVP 0.019 def2-TZVP 0.014 N9 –C9 LANL2DZ 0.013 def2-SVP 0.007 def2-TZVP 0.010 MUE LANL2DZ 0.049 def2-SVP 0.020 def2-TZVP 0.021 a Deviations are computed respect

B97D

ωB97X

ωB97X-D

0.026 0.009 0.010 0.036 0.036 0.036 0.028 0.038 0.040 0.034 0.027 0.034 0.023 0.009 0.020 0.021 0.005 0.016 0.066 0.083 0.067 0.002 0.002 0.003 0.003 0.001 0.001 0.188 0.118 0.123 0.067 0.002 0.010 0.058 0.008 0.004 0.023 0.005 0.015 0.009 0.018 0.002 0.003 0.024 0.008 0.031 0.005 0.008 0.011 0.012 0.009 0.006 0.016 0.014 0.061 0.041 0.043 0.025 0.013 0.013 0.020 0.015 0.014 to the experimental X-Ray data

M06 0.015 0.034 0.039 0.028 0.011 0.006 0.069 0.003 0.005 0.117 0.016 0.008 0.014 0.008 0.022 0.014 0.004 0.010 0.043 0.013 0.015 from Ref.

M06-L 0.043 0.006 0.013 0.024 0.004 0.000 0.058 0.019 0.021 0.116 0.003 0.010 0.007 0.047 0.028 0.022 0.002 0.003 0.045 0.014 0.013 23

N9 –Ru–Cl and C1 –Ru–Cl angles (see Figure 1 for atomic numbering). These geometrical parameters allow us to determine the cymene ring position with respect to the Ru center (through the so-called η 6 binding mode), and to evaluate the covalent bonds established between the metal and the bidentate acridine group. In Tables 1 and 2, we check first the performance of the six functionals employed to optimize the geometry of Complex-I by comparing both the unsigned errors and the mean unsigned errors (MUEs) of the selected bond parameters with respect to the measured X-ray structure. We should note that such a comparison needs to be interpreted cautiously since deviations are not necessarily an indicator of a DFT failure and can be also partially due to the lack of packing effects in the isolated molecule. Nevertheless, the differences between the experimentally recorded and the theoretical predicted geometrical parameters offer a first

10


Page 11 of 39

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


DFT sanity test. As observed in Table 1, the largest errors for the bond lengths obtained using the different density functionals combined to the LANL2DZ basis set, correspond to the Ru–C21 and Ru–C24 bonds for all the functionals, that is, to the two distances which position the cymene ring. The error is especially pronounced for the Ru–C24 bond, ranging from 0.123 ˚ A (ωB97X-D) to 0.116 ˚ A (M06-L). The overall effect in the geometry is partially mitigated by the closer match between crystal data and calculated values for both the Ru– Cl and the Ru–acridine distances (Ru–C1 and Ru–N9 ), eventually resulting in MUE values below 0.061 ˚ A. However, DFT/LANL2DZ calculations cannot be considered accurate enough because the deviation for cymene ring clearly overpasses the acceptable error limit for a Ru– C bond, which should be under 0.07 ˚ A according to the threshold proposed by Calhorda et al. 39 To better handle the most problematic η 6 -cymene bonding mode, we then move on to the def2-SVP basis set. For all functionals, we observe in this case a significant improvement in the positioning of the ring. The errors for the predicted Ru–C21 bond reduce to 0.002–0.003 ˚ A (with the exception of M06-L, 0.019 ˚ A), and the error for the Ru–C24 bond grows slightly, A). Despite the higher A (with the exception of B97D, 0.067 ˚ lying between 0.016 and 0.010 ˚ values provided by the M06-L and B97D functionals, all errors now satisfy the bond distance A). A), resulting in the reduction of the low mean error (MUEs < 0.025 ˚ requirement (

Rh(III) and Ir(III)

Recommend Documents