Machine Learning Driven Optimization of Fe-Based TMCs for Photodynamic Therapy

1 Introduction

Solar energy is a promising alternative to fossil fuels due to its abundance and renewable nature [Kashyap and Ghasemi, 2020], and with the help of photovoltaics it can be converted to electricity [Nazeeruddin et al., 2001] and chemical energy. Transition metal complexes (TMSc), especially those based on nobel metals, are widely used for solar cells, artificial photosynthesis, bioimaging, and photodynamic therapy due to their remarkable photophysical properties [Balzani et al., 2007]. These precious metals, however, are very rare [Xia and Yang, 2017], which stimulates the search for alternative solutions, such as iron. Despite its abundance and lower toxicity, iron is inferior in photophysical properties, because ironbased complexes have lower ligand field and unstable excited states [Wenger, 2019]. Nonetheless, recent progress in prolonging excited state lifetimes in iron-based TMCs opens avenues for efficient photoactive complexes [Kjær et al., 2019].

The photophysics of transition metal complexes plays a key role in their applications. Metal-organic complexes containing π -conjugated ligands have 4 main electronic excited states: ligand-centered (LC), metal-centered (MC), charge transfers from ligand to metal (LMCT), and from metal to ligand (MLCT) [Franc´es-Monerris et al., 2019]. A strong ligand field in Ru complexes results in the destabilization of triplet and quintet MC states, which results in long-lived MLCT states. In the case of iron-containing complexes, the situation is different: due to the weak ligand field, the MC states are lower in energy, which leads to rapid (subpicosecond) deactivation of MLCT and the impossibility of effective use of such compounds in photoactive processes. Recent advances in ligand design, such as the use of strong σ -donor and π -acceptor properties, have extended MLCT lifetimes, making Fe-based complexes viable for photoactive applications. On the other hand, highly symmetric coordination spheres (close to the ideal octahedral disposition) achieved with proper chelating ligands could also destabilize excitedstate deactivation through MC states. We rely on the latter strategy and consider two heteroleptic iron complexes (C1 and C2) (Figure 1) synthetized and characterized by Carrillo and colleagues [Arcos, 2023].

One of the applications of photoactive molecules is PDT, a targeted cancer treatment, which is based on the combined action of a drug and molecular oxygen [Monro et al., 2018]. Upon light irradiation, the excited state can be deactivated through electron or energy transfer to molecular oxygen, which results in reactive oxygen species (ROS). The highly active molecules are reactive towards DNA and cause cell apoptosis. TMCs are very suitable photosensitizers because of their stability, biocompatibility, high absorption within the relative rate of tissue transparency (700-1100 nm) and high triplet state quantum yield [Monro et al., 2018].

Finding the optimal Fe-based photosensitizer for PDT [Bouch´e et al., 2020] can be challenging due to the multiple objectives of the optimization problem and the need for experimental validation. We propose a conceptual ML-driven optimization method to iteratively adjust the molecule. The goal is to achieve a shift in the absorption spectrum toward visible range, thereby enhancing performance in photodynamic therapy, while maintaining key attributes, such as DNA affinity.

2 Proposed Framework

In this section, we present our conceptual ML-driven framework for the iterative optimization of iron-based TMCs, aimed at shifting the absorption spectrum toward the visible range while preserving DNA affinity. While the diffusionmodel-based agent has not yet been implemented, we outline a framework that leverages computationally accessible properties. The framework comprises two core modules: an agent module based on a conditioned diffusion model, and a scoring module, inspired by REINVENT-Transformer [Xu et al., 2023]. The process begins with an initial Fe-based complex (e.g., C1), whose absorption wavelength and DNA binding properties are evaluated using established methods (TDDFT, MD). At the heart of the framework is an agent built on a pretrained diffusion model, adapted from LigandDiff [Jin and Merz Jr, 2024]. The candidate structures are evaluated using TDDFT and MD, and a score is computed as the difference between the desired and the calculated properties. This score guides the agent in generating further modifications. To illustrate the framework, we present results from two ironbased complexes, C1 and C2. C2 is featuring a refined ligand and exhibits a shifted absorption wavelength consistent with our optimization targets, while retaining stable DNA affinity, serving as an example of a model output.

3 Results

3.1 QuantumChemistry

A starting complex is iteratively optimized to shift the absorption wavelength while keeping key PDT features stable. We briefly describe results for complexes C1 and C2.

One of the challenges computational chemists face when studying Fe-based TMCs is the correct determination of the ground state (GS) [Lawson Daku et al., 2005]. For these 100-atom complexes, DFT is suitable, although functionals inconsistently predict the GS. To alleviate these deficiencies, Reiher and coworkers introduced a reparameterization of the B3LYP functional specifically for Fe complexes - B3LYP* functional with reduced (from 20 % to 15%) admixture of the Hartree-Fock exchange [Reiher, 2002; Reiher et al., 2001]. B3LYP and B3LYP* functionals in combination with different basis sets reveal S0 and D1 as the GS for both C1 and C2 systems, respectively. The GAUSSIAN 16 [Frisch et al., 2016] software package has been used in all quantummechanical calculations.

We benchmarked various exchange-correlation functionals and basis sets for C1's key spectroscopic characteristics against B3LYP*. We have studied the following functionals in the benchmark: BLYP, HCTH, M06, M06-2X, PBE0, B3LYP, and CAM-B3LYP, and compiled vertical absorption energies (VAEs), oscillator strengths (f ) and the expectation value of the total spin operator. In the case of the C2 complex, the brightest state D5 lies in the infra-red region, f = 0.1484 and λ = 745 nm (Table 1). The expectation value for the total spin operator is 1.14, which is comparable with an ideal value of 0.75 for a doublet. The results for the C1 complex are similar: the brightest state is again D5, f =0.1614 and λ = 757 nm. The total spin operator is 1.325, indicative of a slightly larger mixing with higher spin states as compared to the bright state of C2.

To understand the nature of the involved excited states, their density matrices obtained with the TD-B3LYP*/631G* method were analyzed utilizing the TheoDORE package [Plasser, 2020]. Results for C2 are presented in Figure 2. The excitation energies and oscillator strengths are given as colored shading. The representation shows that the brightest state is D5, which has a strongly pronounced LMCT character Energy (eV)

Character

3.00

2.00

1.00

0.80

0.60

0.40

0.20

0.00

D2 D5

3.00

3.00

0.00

State D2 D5

D10

State D15

D20

D25

D20

D2 D5 D10 D15 Figure 2: Fragment-based analysis of excited-states in C2 using TheoDORE reveals key LMCT contributions in C2. Molec-

State

State

ular structure (a) of the C2 molecule and its decomposition into

three fragments; (b) excitation energies and oscillator strengths and

excited-state character, automatically assigned using the CT num-

bers

and mixing with LC and LLCT, as shown by the CT numbers. The impact of the pyridine substitution from C2 to C1 on the absorption properties is small.

3.2 Molecular Dynamics

MD simulations rely on a force field to compute forces and equations of motion to describe particle dynamics. The force field must be parametrized, typically against experimental or high-level computational data. We assigned standard AMBER atom types to all ligands and unique atom types to each Fe-ligating atom, and generated the corresponding parameters with VFFDT [Zheng et al., 2016]. All-atom classical MD simulations of complex C2 with a DNA oligonucleotide were conducted to investigate non-covalent interactions by using Amber software [Case et al., 2023]. We visually inspected equilibrium dynamics and used Umbrella Sampling to determine the free binding energy. A 0.5 microsecond equilibruim dynamics simulation showed that C2 interacts with the DNA minor groove for 451.2 ns and remains unbound for 48.8 ns. Figure 3a shows the tert-butyl ligand in the minor groove, with a methyl group oriented toward the DNA. Figure 3b displays the Fe-P distances. Electrostatic forces dominate, as C2's di-cation interacts with the negatively charged phosphate backbone. Small Fe-P fluctuations indicate a stable interaction.

Binding free energy (∆G) was calculated by comparing D20

the free energies of the bound state (15-20 ˚A, where free energy is minimal) and the unbound state (50 ˚A, with no contact). The frequency diagram for 42 windows showed excellent sample overlap, validating the simulation (bias: 1 kcal ·mol-1 · ˚A-2 ). Figure 4's PMF indicates a favorable interaction with an estimated (∆G) of 2 kcal/mol, corresponding to 96 % binding at 310.15 K

The unbound-to-bound transition has a ∼ 1.5kcal/mol barrier (peaking at ∼ 40 ˚A), while the reverse path has a ∼ 3.5 kcal/mol barrier; both are accessible at 310.15 K.

Despite modest binding affinity, C2 remains close to DNA, suggesting that efficiently generated short-lived ROS could reach DNA for PDT.

4 Conclusion

In summary, this study presents an conceptual computational framework that combines quantum chemical calculations, molecular dynamics simulations, and a machine learningdriven optimization module to design and refine Fe-based transition metal complexes for PDT. Our analysis reveals that the near-IR absorption is primarily governed by the LMCT states. Furthermore, the molecular dynamics simulations indicate that while the complexes exhibit modest binding affinities to the DNA, they remain in close proximity, a key feature for effective ROS mediated damage in PDT.

The proposed ML-driven framework, though not yet implemented, demonstrates a feasible approach to iteratively tailor ligand structures by leveraging rich simulation data. Our work establishes a foundation for future optimization, showing that key PDT-relevant properties can be computationally evaluated and used to guide molecular design.

Overall, our work not only advances the fundamental understanding of the photophysical and binding properties of Fe-based complexes but also lays the foundation for future studies aimed at optimizing these compounds for practical PDT applications.

Acknowledgements

The work of Ilya Makarov on Section 1 was supported by the Ministry of Economic Development of the RF (code 25-13966879-1-0003).

0.00