{'fraction_non_alphanumeric': 0.03648692282854375, 'fraction_numerical': 0.021059806981666845, 'mean_word_length': 4.999569522169608, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Singlet Oxygen Generation by Porphyrins and Metalloporphyrins Revisited: a Quantitative Structure-Property Relationship (QSPR) Study

Abstract Porphyrins and metalloporphyrins are used as photosensitizers in photocatalysis, photodynamic therapy (PDT), disinfection, degradation of persistent pollutants and other applications. Their mechanism of action involves intersystem crossing to triplet excited state followed by formation of singlet oxygen (1O2), which is a highly reactive species and mediates various oxidative processes. The design of advanced sensitizers based on porphyrin compounds have attracted significant attention in recent years. However, it is still difficult to predict the efficiency of singlet oxygen generation for a given structure. Our goal was to develop a quantitative structure-property relationship (QSPR) model for the fast virtual screening and prediction of singlet oxygen quantum yields for pophyrins and metalloporphyrins. We performed QSPR analysis of a dataset containing 32 compounds, including various porphyrins and their analogues (chlorins and bacteriochlorins). Quantum-chemical descriptors were calculated using Density Functional Theory (DFT), namely B3LYP and M062X functionals. Three different machine learning methods were used to develop QSPR models: random forest regression (RFR), support vector regression (SVR), and multiple linear regression (MLR). The optimal QSPR model «structure – singlet oxygen generation quantum yield» obtained using RFR method demonstrated high determination coefficient for the training set (R2 = 0.949) and the highest predicting ability for the test set (pred_R2 = 0.875). This proves that the developed QSPR method is realiable and can be directly applied in the studies of singlet oxygen generation both for free base porphyrins and their metal complexes. We believe that QSPR approach developed in this study can be useful for the search of new poprhyrin photosensitizers with enhanced singlet oxygen generation ability.


Introduction
Organic photosensitizers, efficiently forming long-living triplet excited states and singlet oxygen ( 1 O2), found applications in various fields of technology, e.g. photocatalysis, disinfection, photodynamic therapy (PDT) and degradation of persistent pollutants [1]. In recent years there has been a growing number of studies devoted to synthesis and photophysical investigation of novel photosensitizers [2]. However, despite a great number of papers on this topic it is still challenging to design a photosensitizer with predetermined photophysical parameters, in particular absorption in a target region of spectrum, high triplet excited state yield and long lifetime. For this reason, development of structure-activity relationships for predicting photosensitizing ability of organic molecules is crucially important.
Porphyrins represent a unique class of organic photosensitizers due to a combination of easily controllable photophysical properties and synthetic accessibility [3]. Porphyrins have multiple applications including drug delivery [4,5], bioimaging [6,7], and non-linear optics [8,9] to name a few. One of the most intriguing properties of porphyrins is their ability to produce singlet oxygen with high quantum yield under visible light irradiation, which is used for the treatment of skin diseases such as psoriasis [10], port-wine stains [11], as well as different types of cancer like mycosis fungoides (a type of a cutaneous T-cell lymphoma) [12], malignant mesothelioma [13], glioma [14], brain tumor [15], fibrosarcoma [16], colorectal cancer [17], osteosarcoma [18].
Porphyrins and metalloporphyrins are used in antimicrobial PDT [19][20][21]. Moreover, some of the photosensitizers are used for the photodynamic treatment of such skin diseases as vitiligo, psoriasis, etc.
In cells, porphyrins may participate in photooxidation reactions of electron donors (nucleic acids, peptides, lipids, and other biomolecules). This type of reactions relate to type I photooxidation mechanism: P (S0) → 1 P* → 3 P* (eq. 1) As one can see both types of photooxidation reactions depend on the triplet state formation. We may see that an efficient porphyrin photoagent possesses: 1) good penetration into the cells; 2) high absorption coefficient in the red part of the visible spectrum; 3) high triplet state quantum yield; 4) long triplet state lifetimes; 5) high ROS generation efficiency, in particular, high quantum yield of singlet oxygen generation (Φ∆).
The parent compound of porphyrins, porphine (structure I, Fig. 1), exhibits efficient intersystem crossing (ISC). Chemical modification of porphine scaffold allows to optimize its properties for practical applications. For instance, chlorin (II) and bacteriochlorin (III) derivatives possess intense absorption in the near-IR window of tissue and high values of Φ Δ , that makes them ideal photsensitizers [22]. Introduction of peripheral substituents into meso and  positions of the macrocycle can be used to modulate the wavelength of absorption maxima within the therapeutic window and tune the solubility. However, it is hard to predict the effect of substituents on ISC and in some cases functionalization (e.g. with strong electron accepting groups) can lead to reduced triplet state yields and lifetimes [23]. Quantitive Structure-Activity Relationship (QSAR) models are commonly used for the prediction of absorption, distribution, metabolism, excretion (ADME), toxicity, and selectivity towards certain cells in pharmacokinetics. To our knowledge QSAR models for the prediction of phototoxicity and for 1 O2 generation quantum yield prediction are rare. However, earlier we showed that Quantitive Structure-Property Relationship (QSPR) method may be used for the prediction of singlet oxygen generation quantum yield for organic sensitizers [24,25].
The quantum yield of singlet oxygen generation Φ∆ is a ratio between the amount of generated 1 O2 molecules and the quantity of photons absorbed. Several parameters influence the generation of 1 O2: steric and structural effects, triplet state energy (ET), oxidation potential, electronic configuration, solvent polarity, aggregation, and excimer formation [26]. Interaction both with singlet and triplet excited states of a sensitizer can lead to energy transfer. However, for most photosensitizers, the role of ground state molecular oxygen interaction with S1 state of the photosensitizer is limited to the enhancement of the triplet quantum yield, and the main pathway of 1 O2 generation is O2 interaction with the first excited triplet state, T1 [26]. Thus, Φ∆ is determined as: where ΦT is the triplet state formation quantum yield, pT O2 is the fraction of T1 state quenched by molecular oxygen, fT Δ is the efficiency of 1 O2 production during molecular oxygen quenching of T1. Because the intramolecular T1 to S0 transition is spin-forbidden, T1 states generally have long lifetimes, allowing almost complete quenching by O2 (pT O2 ≈ 1) in air-equilibrated solutions. Incomplete quenching is observed only with triplet states possessing a short lifetime due to high chemical reactivity or fast intersystem crossing (ISC) (i.e., nπ* triplets), or in the presence of very high concentrations of other quenchers [26].
The higher the rate of physical and chemical quenching of the triplets, the lower the efficiency of energy transfer from the triplet state sensitizer to O2 [27]: where ket is the rate of energy transfer reaction, kq T is the rate constant of other possible processes involving 3 P* and O2 (e.g. superoxide anion formation, quenching of 3 P* by molecular oxygen leading to the ground state of both reactants); kp and kISC are the rate constant of radiative decay of the triplet state photosensitizer 3 P* and the rate constant of non-radiative decay, respectively. It was found that magnitudes of kq T and f T Δ are inversely correlated, and both parameters exhibit pronounced sensitivity to the oxidation potential [28].
The quantum yield of the triplet state formation depends on the rate of intersystem crossing k'ISC and the rates of competing processes: where k'ISC is the rate constant of S-T intersystem crossing, kq is the rate constant of physical and chemical quenching of singlet excited state, kfl is the fluorescence rate constant, kic is the internal conversion rate constant (Fig. 2). ISC is governed by spin-orbit coupling (SOC) [29], which is a relativistic effect. SOC is enhanced in the presence of atoms with large nuclei (heavy atoms), such as transition metals and halogens (bromine and iodine). Introduction of heavy atoms into organic chromophores is a common strategy used to access efficient photosensitizers. However, incorporation of heavy atoms can lead to other problems, such as higher costs, toxicity, lower photostability and also a shortened triplet state lifetime [30]. These issues have stimulated much efforts in the field and photosensitizers with efficient ISC in the absence of any heavy atoms are being actively investigated [31,32].
There are qualitative rules, describing ISC. According to El-Sayed rule [33], S-T ISC is efficient between 1 nπ* and 3 ππ*, 1 ππ* and 3 nπ* states. Another generally known rule is the energy gap law which says that transition probability depends exponentially on the adiabatic energy difference ΔE, i.e., the smaller the energy gap the larger the transition probability [34]. In quantum chemistry approximate SOC operators and stateof-the-art theoretical methods are used for evaluation of ISC rates. Since for pharmaceutical purposes it is required to perform virtual screening on large libraries of compounds, calculation of SOC operators for these purposes seems inappropriate. The alternative strategy proposed in this study is the application of the QSPR and machine learning methodology for the prediction of Φ Δ . Photosensitizer synthesis and experimental determination of singlet oxygen generation quantum yield is a time-and labor-consuming process. Therefore, preliminary, pre-synthetic screening of sensitizers for their ability to generate 1 O2 would be of great value. That is why our goal was to build QSPR models predicting the quantum yield of 1 O2 generation.
To develop a reliable QSPR model, we focused on a series of compounds, combining several classes of functionalized porphyrins and their analogues (Fig. 1b). It includes four sub-groups: 1) tetraphenyl porphyrins derivatives with different substituents in meso-phenyl rings and -positions of pyrrolic rings; 2) -alkyl substituted porphyrins; 3) derivatives of etioporphyrin and deuteroporphyrin bearing electron donating or accepting groups; 4) chlorin and bacteriochlorin derivatives. In addition, zinc (II) complexes of six structurally different porphyrins were included into the dataset.
Previously reported data on ΦΔ of these compounds in organic solvent were used for analysis. We performed QSPR analysis and optimized the model for predicting photosensitizing ability. We believe that the results of our study will be helpful for the virtual screening of new synthetic porphyrins and metalloporphyrins.

Methods
All porphyrin compounds studied in this work are poorly soluble in water, but possess good solubility in organic solvents. For this reason, experimental values of ΦΔ measured in toluene were used for analysis [35]. The activity of compounds, presented as Log ΦΔ, was used as a dependent variable in the QSPR study. The porphyrins were divided into the training (80% of studied compounds) and test sets (20%) using random number generation. The requirements for the maximum and minimum values of Log ΦΔ in the test set were the following: 1) the maximum Log ΦΔ value should be less than or equal to the maximum value of the training set; 2) the minimum Log ΦΔ value should be higher than or equal to the minimum value of the training set.
Two datasets were produced from the original one: 1) to simplify the task 26 free base porphyrins were included in the first dataset; 2) the second one consisted of free base porphyrins and corresponding zinc (II) porphyrins (32 molecules in total). Chemical structures of studied compounds are presented in the Supporting Information (Fig. S1).

Geometry optimization
Conformational analysis of porphyrin molecules was performed using Spartan v. 16 modeling software from Wavefunction, Inc. (www.wavefun.com). Generation of low-energy conformers was performed using PM6 semi-empirical quantum-chemical method [36].
Geometry optimization was done using Density Functional Theory. We used two different methods for geometry optimization and molecular descriptors calculations. B3LYP [37,38] functional with 6-31G(d,p) basis set for all atoms and LANL2DZ basis set with effective core potential (ECP) [39][40][41] was used for iodine in Spartan v. 16. We used M06-2X [42,43] with mixed basis set of 6-31G(d,p) for all atoms and SDD basis set with ECP for iodine [44] in Gaussian v. 16 [45]. Optimized porphyrin metalloporphyrin geometries are presented on Fig. S2.

Molecular descriptors
Dragon v. 7.0 descriptors were calculated using online chemical modeling environment

Machine learning techniques. Model search
Machine learning (ML) was performed using scikit-learn library of Python programming language. We used three different ML methods: support vector regression (SVR), multiple linear regression (MLR), and random forest regression (RFR) (Fig. 3). We performed standard-scaling of descriptor values for SVR models. The idea of standard-scaling is to make the values of each descriptor in the dataset have zero mean and unitvariance, according to the expression:

Descriptors contribution to the model
We estimated the contribution of each descriptor to the MLR models using following equation: where α(x1) is the relative contribution of the descriptor x1 to the model with several descriptors, 2 ( ) is the determination coefficient of the nth descriptor towards Log ΦΔ.
Relative contribution of descriptors in SVR and RFR models was evaluated by using the facilities implemented in these methods: feature importance was evaluated in RFR; in SVR with linear kernel coefficients of the descriptors were used for relative contribution calculation according to the following formula: where ( 1 ) is the value of the coefficient for descriptor x1; ∑ (" + ") and ∑ (" − ") are the sums of coefficients for all positive and negative coefficients, respectively.

Statistical parameters
Best QSPR models were selected on the basis of statistical parameters such as R 2 (coefficient of determination of the training set), q 2 (internally cross-validated R 2 using the leave-one-out (LOO) method), and pred_R 2 (predictive R 2 for the test set of compounds). Q 2 and pred_R 2 were used for model validation and comparison. QSPR models obtained were validated and tested for predicting ability using a test set of compounds.
To calculate q 2 parameter each molecule in the training set was excluded once and Log ΦΔ of the excluded molecule was predicted by using the model developed by the remaining compounds. Q 2 describes the internal stability of a model and was calculated using the following formula: where yi and ŷi are the actual and predicted Log ФΔ value of the ith molecule in the training set, respectively; ymean is the average Log ФΔ of all compounds in the training set.
For external validation, the activity of each compound in the test set was predicted using the model developed by the training set. Pred_R 2 indicates the predictive ability of the model.
Pred_R 2 was calculated according to the equation: where yact and ypred are the actual and predicted activity of the ith compound in the test set, respectively; ymean is the average Log ФΔ of all molecules in the training set. Both summations are over all compounds in the test set.

Results and Discussion
To simplify the task we excluded metalloporphyrins and started with the analysis of a dataset consisting of 26 free base porphyrins [35]. Most of the compounds in this dataset are derivatives of tetraphenylporphyrin (TPP), containing different substituents in meso-phenyl groups. We also included various representatives of -alkyl-substituted porphyrins as well as derivatives of chlorin and bacteriochlorin. We used three different machine learning methods in our analysis: random forest (RFR), support vector regression (SVR), and multiple linear regression (MLR). Statistical parameters of the models obtained are summarized in Table 1.

Model 1 (support vector regression)
Let's regard the models obtained. Thus, we built a SVR model with linear kernel (Model 1). The model possesses constants C and epsilon both equal to 0.0025. Linear kernel makes it possible to evaluate descriptors contribution. The SVR model uses a full set of 182 descriptors, but the two most influential descriptors are MATS1m (relative contribution -1.7%) and HATS4m (1.5%). MATS1m is a Moran autocorrelation lag 1 / weighted by atomic masses (2D autocorrelation descriptor) [48]. The descriptor is inversely correlated with Log ФΔ. N-Methyloctaethylporphyrin has the highest value of MATS1m equal to 0.063, while tetra(4bromophenyl)porphyrin ((p-Br)4-TPP) has the lowest value equal to -0.002. In general, halogenated molecules have low values of the descriptor while halogen-free molecules have high values of MATS1m. That means that halogenation is favorable for high ΦΔ, which is a wellknown effect [49]. HATS4m is a leverage-weighted autocorrelation of lag 4 / weighted by atomic masses (GETAWAY descriptor) [50,51]. The descriptor is directly proportional to Log ΦΔ. The value of HATS3m is equal to zero for 10 compounds presented in the dataset (   Table S1 and in Fig. 4. We randomly divided the dataset into the training and testing sets four times and tested the models to additionally confirm good quality of the model obtained. Thus, key parameters had the following value ranges after additional dataset splitting: R 2 = 0.952-0.993, pred_R 2 = 0.610-0.691, q 2 = 0.524-0.734, and RMSE = 0.019-0.037. These results show that Model 1 is tolerant to alternative dataset splitting, possesses high predicting ability and internal stability. MATS1e (relative contribution -43.4%) is a Moran autocorrelation lag 1 / weighted by atomic Sanderson electronegativities [48]. MATS1e is inversely correlated with Log ФΔ. Among the studied compounds 5,10,15-NO2-etioporphyrin I has the highest value of MATS1e equal to 0.279, while porphine has the value of the descriptor equal to -0.086 (Table 3). Electronegative nitrogen atoms are evenly distributed in porphine, which makes its MATS1e value low, while in 5,10,15-NO2-etioporphyrin I, 5,10-NO2-etioporphyrin I, and 5,15-NO2-etioporphyrin I introduction of nitro groups results in asymmetrical distribution. Apparently, the main point is that MATS1e descriptor allows to distinguish molecules bearing electron accepting groups, e.g. nitro groups, from molecules lacking such groups. Tetraphenyltetrahydroporphyrin has the highest value of S1_fOSC equal to 0.4625 while (o-Cl)4-TPP has the lowest value of the descriptor (0.015). According to Kasha's rule, emission always occurs from the lowest excited state of the same multiplicity (S1). The oscillator strength of a transition between the ground state and the excited state is related directly to the intensity of emission [54,55]. Apparently, molecules with higher S1_fOSC possess higher emission rate and lower intersystem crossing rate (a competing process), which results in the decrease of both ФT and ФΔ.
Experimental and predicted by the MLR model ΦΔ values for each porphyrin compound are presented in Table S2 and in Fig. 4.
We randomly splitted the dataset four times and tested the model performance in each case. Key statistical parameters had the following value ranges: R 2 = 0.844-0.922, pred_R 2 = 0.646-0.905, q 2 = 0.495-0.827, and RMSE = 0.034-0.040. As one can see, statistical parameters after alternative dataset splitting are sufficient.

Model 3 (random forest regression)
Next, we have developed a random forest regression model (Model 3). The model selected 54 descriptors out of the total of 182 descriptors; eight trees with a maximal depth equal to 10 were used by the RFR model. Two descriptors have a relative contribution more than 10%.
We regarded these descriptors in details.   Table S3 and in Fig. 4.
After getting a proof-of-principle for QSPR application in the study of singlet oxygen generation by free base porphyrins, we decided to extend this method to a larger dataset of molecules, which includes corresponding metalloporphyrins. A dataset consisting of twenty-six free base porphyrins and six zinc (II) porphyrins was analyzed.
Statistical parameters of the models obtained are summarized in Table 5. As one can see from Table 5, the statistical parameters of all three models are sufficient.
The SVR model has the best values of R, R 2 , q 2 , root mean square error, and root mean square ΦΔ error. The RFR model has the best values of maximum ΦΔ error and pred_R 2 . Pred_R 2 is probably the most important parameter evaluating the predictive ability of the model towards the test set. We suppose that RFR model is the optimal QSPR model «structuresinglet oxygen generation quantum yield» since it possesses high determination coefficient for the training set (R 2 = 0.949) and the highest predicting ability towards the test set (pred_R 2 = 0.875).
Key statistical parameters after alternative splitting had the following value ranges, is barely sufficient (0.495). Nevertheless, we suppose that all three models possess satisfactory tolerance to alternative dataset splitting.

Model 4 (SVR)
We built a support vector regression (SVR) model with linear kernel (Model 4). The model possesses constants C and epsilon both equal to 0.005. Linear kernel allows to evaluate relative contributions of the descriptors in the model. The SVR model uses a full set of 182 descriptors, but the two most influential are S1_energy (relative contribution 3.1%) and Eig02_EA(dm) (-2.5%). S1_energy is an energy of the first electronic transition S0S1 calculated with M06-2X/6-31G(d,p) method for the gas phase. The descriptor is directly proportional to Log ФΔ. Zn-Etioporphyrin possesses the highest value of S1 state energy (2.437 eV) while octamethyl-TPP has the lowest value of S1_energy equal to 2.036 eV. Hence, the importance of S0S1 transition arises again: see Model 2 for comparison. We suppose that the role of S1 state energy reveals a well-known phenomenon that higher S1 energy is proportional to a lower rate of chemical quenching of 1 O2 [26]. the compounds capable of 1 O2 quenching either do not generate singlet oxygen or do it less efficiently. Eig02_EA(dm) is an eigenvalue No 2 from edge adjacency matrix weighted by dipole moment, which was previously used in Model 2 and allows to distinguish molecules with strong electron accepting groups.  Table S4 and in Fig. 5.  Fig. 5. Experimental vs. predicted Log ФΔ for the dataset combining free base porphyrins and metalloporphyrins.
We randomly divided the dataset into the training and testing sets four times and tested the models to additionally confirm good quality of the Model 4. Main statistical parameters had the following value ranges after additional dataset splitting: R 2 = 0.976-0.985, pred_R 2 = 0.567-0.881, q 2 = 0.654-0.799, and RMSE = 0.016-0.031. These results show that Model 4 is tolerant to alternative splitting and possesses high predicting ability and stability: determination coefficients R 2 are especially high.

Model 5 (MLR)
We MATS1e (relative contribution -36.9%) is a Moran autocorrelation lag 1 / weighted by atomic Sanderson electronegativities [48]. MATS1e is inversely correlated with Log ФΔ. As well as Eig02_EA(dm), MATS1e allows to distinguish molecules bearing electron accepting groups.  Table 7). The diagonal elements are roughly proportional to the polarizability of the atoms. The larger the polarizability of the atoms the larger the SpMax4_Bh(p) value, which leads to a increase in Log ФΔ.

Model 6 (RFR)
Next, we have built a random forest regression model (Model 6). The model used 70 descriptors and 12 trees with a maximal depth equal to 80. Two descriptors had a relative contribution more than 10%. We will regard these descriptors in details.  We also calculated vertical S1 and T1 energies of all porphyrin compounds. Surprisingly, the singlet-triplet energy gap was not found to play a role in any of the models. In fact, the determination coefficient R 2 between S-T gap values and Log ΦΔ was equal to 0.02. This is consistent with our previous work on furocoumarins, which demonstrated that S-T gap does not influence the Log ΦT prediction [25].

SpDiam_EA
Key statistical parameters after alternative splitting had the following value ranges, Thus, we proved that QSPR method is applicable to the study of singlet oxygen generation not only for free base porphyrins, but also for metalloporphyrins. The developed models deliver reasonable results and are theoretically well interpretable.

Conclusion
In this study, a QSPR approach for predicting the efficiency of singlet oxygen generation by free base pophyrins and metalloporphyrins was developed. Three different methods of machine learning were applied: support vector regression (SVR), multiple linear regression (MLR), and random forest regression (RFR).
First, a dataset consisting of 26 free base porphyrins was studied and three models predicting ΦΔ values were developed. All three models were found to be sufficient. The SVR model was found to be most accurate in the prediction of ΦΔ for the training set of compounds. The MLR model showed the best performance for predicting ΦΔ for the test set. The RFR model was found to be the most internally stable one (q 2 = 0.787). Thus, we proved that QSPR methodology is applicable to the study of 1 O2 generation by porphyrins.
Second, we examined the activity for a larger set of molecules, comprised of 26 free base porphyrins and six metalloporphyrins. All three models were found to be sufficient. The SVR model performed as the best one in predicting ΦΔ for the training set and was the most internally stable one (q 2 = 0.842). The RFR model showed the best performance in predicting ΦΔ for the test set (pred_R 2 = 0.875). The MLR model showed the best tolerance to alternative dataset splitting. Thus, we showed that QSPR methodology is applicable to the study of 1 O2 generation by coordination compounds, in particular, metalloporphyrins.
The models were built using a combination of autocorrelation, GETAWAY, topological, and quantum-chemical descriptors. SpDiam_EA(bo), Eig02_EA(dm), MATS1e, and S1 state energy descriptors were found to be most influential. Relative contribution of the descriptors to the models was evaluated and interpreted taking into account structural and electronic effects in the molecules. MLR models provided valuable information about the descriptors involved. For example, it proves that halogenation of porhyrins increases ΦΔ, while introduction of nitro groups has a negative influence on ΦΔ.
Machine learning and QSPR are largely used nowadays in photochemistry. For example, they are applied for the prediction of photovoltaic properties of new materials [67,68]. When we began our studies in 2016 [24], there were almost no works dedicated to 1 O2 generation analysis with QSPR, but now our studies belong to a new emerging field of knowledge [69]. The developed models can be used not only for PDT, but also for the PUVA-therapy of such diseases as psoriasis and vitiligo [25]. We believe that the QSPR models developed in our study will be useful for the screening and design of free base porphyrin and metalloporphyrin photosensitizers with high singlet oxygen generation quantum yields.