In this study, all experimental data of logP_oct, logP₁₆ (the logarithm of the partition coefficient between hexadecane and water), logP_chl (the logarithm of the partition coefficient between chloroform and water), logP_aln (the logarithm of the partition coefficient between aniline and water), logK_brain (the logarithm of the partition coefficient from air to human brain) and logK_p (logarithm of experimental human skin permeability) are collected from literatures (Abraham et al., 1994; Abraham et al., 1999; Abraham and Martins, 2004; Abraham et al., 2006; Abraham et al., 2015; Zhang et al., 2017). Hydrogen bond acceptors (HBAs) include very weak H-bond acceptors. For example, the sp2 carbon atoms from carbon-carbon double bonds and aromatic rings are weak HBAs. Hydrogen bond donors (HBDs) include very weak H-bond donors. For example, the hydrogen atoms in CHCl₃ and CH₃NO₂ are weak HBDs.

S_m is a molecular descriptor developed in this study. The S_m values of organic compounds were calculated based on the formula of the compounds. Assume the formula of a neutral organic compound is C_cH_hO_oN_nS_sF_fCl_clBr_brI_i, the S_m of this compound is S m =c+0 .3h+o+n+2s+0 .6*f+1 .8cl+2 .2br+2 .6i-0 .2 N c3 -0 .6 N c4 . (1)

In Eq. 1, c, h, o, n, s, f, cl, br and i are the numbers of carbon, hydrogen, oxygen, nitrogen, sulfur, fluoride, chloride, bromide and iodide atoms in the solute, N_c3 is the number of sp3 carbons connecting three heavy atoms (fluoride is not included), N_c4 is the number of sp3 carbons connecting four heavy atoms (fluoride is not included).

H_{M_HBD} values of solutes were calculated based on the approach reported in a previous study (Chen et al., 2020).

In this study, the flexibility of a solute is calculated by summarizing the flexibilities of the bonds of the solute. If a bond is not rotatable or if the rotation of a bond does not change the conformation of the solute, the flexibility of the bond is set to zero (note: hydrogen atoms are not included for determining conformations). The flexibility of the C—C bond in R¹CH₂—CH₂R² is set to one. If the energy barrier for rotating a bond is obviously higher than that for rotating the R¹CH2—CH₂R (Sun et al., 2019) bond, the flex value is set to zero. For example, the C—N bond in RCO—NH and the C—C bond in Ar—CO are set to zero. If the energy barrier for rotating a bond is obviously lower than that for rotating the R¹CH₂—CH₂R² bond, the flex value is set to 1.5. For example, the energy barrier for rotating the R¹O—CH₂R² bond is lower than that for R¹CH₂—CH₂R (Sun et al., 2019) and thus the Flex value of the C—O bond is set to 1.5. Also, the flexibility of C—C in R¹CH₂—C₆H₅ is set to 0.5 because of the symmetry of phenyl ring.

The free energy changes for transferring depolarized solutes from water to hexadecane (ΔG_{tr_depol}) were calculated based on the method reported in previous study (Chen et al., 2020). Based on the logP_oct (or logP_chl) and ΔG_{tr_depol} values of nonpolar compounds, the model for the regression of logP_oct (or logP_chl) against ΔG_{tr_depol} was developed. This model was then used to calculate the logP_oct (or logP_chl) values for depolarized solutes. For a solute containing HBAs but no HBDs, the difference between the calculated logP_oct (or logP_chl) for the depolarized solute and the experimental logP_oct (or logP_chl) of the solute is the effect of HBAs on the logP_oct (or logP_chl) of the solute.

All the models and the statistical reliabilities of the models were obtained by performing the multiple linear regressions implemented in Excel.

In the theoretical derivation, we used “Y” to represent a property and the symbol “ΔG_Y” to represent the free energy change that determines Y. The ΔG_Y values for many properties are not easy to be calculated directly. Thus, we decomposed ΔG_Y into the free energy changes that are caused by the factors affecting Y. The free energy change caused by a factor is denoted by “ΔG_F”. Thus, ΔG_Y equals the summarization of the ΔG_Fs for all the factors affecting Y. Δ G Y = ∑ Δ G F (2)

For the properties depending on the noncovalent interactions of solutes, they are affected by the molecular sizes, hydrogen-bond acceptors (HBAs), and hydrogen-bond donors (HBDs) of the solutes, which was demonstrated in a previous study (Chen et al., 2020). Many properties are also affected by the flexibilities of solutes. For example, the partition coefficients of organic compounds between a flexible environment (e.g., blood) and a much less flexible environment (e.g., muscle) are obviously affected by the flexibilities of the compounds. It is challenging to accurately quantify the ΔG_Fs for various properties. However, it is possible to develop molecular descriptors that are directly proportional to ΔG_Fs. We used D_F to represent the molecular descriptor that is directly proportional to ΔG_F. Then, ΔG_Y can be expressed as: Δ G Y = ∑ k F D F (3)

The k_F values are constant for a given property. Theoretically, if the molecular descriptors apply to various properties, Eq. 3 can be used to construct models with high predictive power for the properties. Many properties are mainly affected by the molecular sizes, HBAs, HBDs and flexibilities of solutes. Thus, in this study, we developed or selected molecular descriptors for quantifying the effects of molecular size, HBAs, HBDs and flexibility on the properties.

The molecular descriptor we developed for quantifying the effects of molecular size on properties is denoted by “S_m”. The S_m values of organic compounds represent the relative molecular sizes of the compounds and can be easily calculated from their molecular formulas (see Computational Methods). For example, the S_m for catechol (formula: C₆H₆O₂) is 9.8 (num. for C + 0.3 ×   num. for H + num. for O). To illustrate whether S_m is an ideal molecular descriptor for molecular size, we first explored the linear associations between logP₁₆ and S_m and between logP_oct and S_m for a series of alkane compounds (Figure 1A). The logP₁₆ and logP_oct values for alkane compounds are affected merely by the sizes of the compounds. The robust linear associations in Figure 1A support that S_m is directly proportional to the effects of molecular size on logP₁₆ and logP_oct. We next explore whether S_m is also an ideal molecular descriptor of molecular size for the properties that have little correlation with logP₁₆ or logP_oct. As reported in a previous study, logK_brain has little correlation with logP₁₆ and logP_oct (Chen et al., 2020). We thus explored the linear association between logK_brain and S_m for nonpolar solutes (Figure 1B). The R² and SD values indicate that there is a strong linear association between logK_brain and S_m. In Figure 1C, we plotted the free energy changes for transferring the depolarized compounds from water to hexadecane (ΔG_{tr_depol}) against the S_m values for the compounds from Supplementary Table S1 of a previous study (Chen et al., 2020). The high statistical reliability for the regression of ΔG_{tr_depol} against S_m further supports that S_m is an ideal molecular descriptor for quantifying the effect of molecular size on the properties depending on noncovalent interactions. Thus, S_m is an ideal molecular descriptor for molecular size and applies to various properties.

In previous studies, we defined the water to hexadecane phase transferring free energy contributed by the electrostatic interactions of the HBAs of a solute as the overall H-bonding capability of the HBAs of the solute (Chen et al., 2019; Chen et al., 2020; Chen et al., 2016) and this overall H-bonding capability is donated by “H_{M_HBA}.” The definition indicates that H_{M_HBA} is an ideal molecular descriptor for quantifying the effects of HBAs on logP₁₆. We next explored whether H_{M_HBA} is an ideal molecular descriptor for logP_oct and logP_chl. The strong linear associations between the effect of HBAs on logP_oct and H_{M_HBA} (Figure 2A) and between the effect of HBAs on logP_chl and H_{M_HBA} (Figure 2B) suggest that H_{M_HBA} is an ideal molecular descriptor for quantifying the effects of HBAs on various properties. Similarly, we defined the water to hexadecane phase transferring free energy contributed by the electrostatic interactions of the HBDs of a solute as the overall H-bonding capability of the HBDs of the solute (H_{M_HBD}) (Chen et al., 2016; Chen et al., 2019; Chen et al., 2020). In a previous study, we revealed that the contribution of a protein-ligand H-bond to the protein-ligand binding free energy is directly proportional to the H-bonding capability of the HBA and the H-bonding capability of the HBD (Chen et al., 2016). We also found that the effect of an enzyme-substrate H-bond interaction on the free energy barrier of the enzymatic reaction is directly proportional to the H-bonding capability of the atom from the enzyme (Chen et al., 2019). Thus, we believe that the effects of HBAs and HBDs of solutes on the properties related to noncovalent interactions are directly proportional to the H_{M_HBA} and H_{M_HBD} values of the solutes. H_{M_HBA} and H_{M_HBD} are ideal molecular descriptors for quantifying the effect of HBAs and HBDs on the properties related to noncovalent interactions.

The molecular descriptor for quantifying the effect of molecular flexibility on properties is denoted by “Flex.” The effects of molecular flexibility on properties mainly result from rotatable bonds of the solutes because the rotatable bonds of the solutes can rotate more freely in some environments than in other environments. The flexibilities of solutes are calculated from the rotatable bonds of the solutes, especially the rotatable bonds that change the conformations of solutes (see Computational Methods). Thus, for many properties that are affected by molecular size, HBAs, HBDs and flexibility, they can be quantified with the following equation Y= k 1 S m + k 2   H M_HBA + k 3 H M_HBD + k 4 Flex+c (4) where k₁, k₂, k₃, k₄ and c are constants for a give property. Organic compounds usually contain multiple HBAs and the HBAs affect each other. The accurate calculation of H_{M_HBA} for many organic compounds is not easy. Eq. 4 would become simpler and easier to use if H_{M_HBA} is replaced by logP_ow because logP_ow is a well-known molecular descriptor for predicting properties (Liu et al., 2019; Zhu et al., 2020b) and can be obtained accurately via experimental and/or computational approaches. Based on the fact that logP_ow is a property and Eq. 4 also applies to logP_ow, we can convert Eqs. 4 to 5 (see Supplementary Text S1 for the detail of the process of the conversion). Y= b 1 log P ow + b 2 S m + b 3 H M_HBD + b 4 Flex+c (5) where b₁, b₂, b₃, and b₄ are constants, logP_ow is logP₁₆ or logP_oct. Eq. 5 is identical to Eq. 4. Both equations are correct for the properties that are determined by the noncovalent interactions of solutes with flexible environments. All the factors related to effects of noncovalent interactions on phase-transferring free energies, including electrostatic interaction, desovation, van der Waals interactions, entropy change, etc. are considered in Eq. 5. Eq. 5 is the general LFER we developed for predicting the properties that depends on the noncovalent interactions of solutes with flexible environments. Although S_m and H_{M_HBD} may be strongly correlated with logP_ow for some properties, none of the molecular descriptors can be omitted because Eq. 4 is a general LFER for various different properties.

Because logP_oct and logP₁₆ can be calculated accurately from the structures of solutes (Chen et al., 2020), it is expected that this method still performs well when all of the molecular descriptors in the LFER are calculated from solute structures. For example, the R², Q² _ext and SD of the model for predicting human skin permeability, in which all predictive variables are calculated from solute structures, are 0.940, 0.957, and 0.202 (see Supplementary Text S4). Thus, the general LFER developed in this study has obvious advantages in predicting many properties related to noncovalent interactions.

Above examples indicate that the models constructed according to the LFER for many specific properties have high predictive power. Moreover, the performance of the models is independent of the compounds for investigation, suggesting that the models can provide guidance for improving properties of organic compounds and designing compounds with optimal properties. The merits of the LFER result from the theoretical derivation, which ensures that the quantitative relationships in the models constructed according to the LFER are correct in the aspect of thermodynamics. For the QSPR models developed using mathematical and statistical tools, the predictive variables are selected from a few thousand molecular descriptors based on the data of the compounds in training sets. The relationships between the properties and molecular descriptors in the QSPR models are statistical relationships for the compounds in training sets. The QSPR models usually work well only for the compounds in the training set and similar compounds, but may do not work well for other compounds. Thus, for the properties determined by the noncovalent interactions of solutes with flexible environments, the models developed according to the proposed LFER performs better than the QSPR models developed by using mathematical and statistical tools, including robust artificial neural networks. Developing models according to the proposed LFER is faster and computationally cheap than developing traditional QSPR models because the process of the variable selection is not required. Moreover, the proposed LFER is quite simple and can be easily used by the researchers across various fields, while expert knowledge is required for developing robust artificial neural networks, such as the knowledge in choosing the most appropriate approach. Thus, the method developed in study has obvious advantages over the traditional QSPR construction method. Thermodynamics-based theoretical derivation can be used to solve many problems that are hard to be solved by using mathematical and statistical tools. In addition, results in this study demonstrate that there are quantitative relationships between the properties related to thermodynamics, suggesting that many properties can be accurately predicted from other properties.

The theoretical derivation in this study is based on the assumption that solutes have similar interactions with their environments, which requires that the environments for the properties are flexible or the properties have little relationship with the conformation or orientations of solutes. Thus, the present LFER may not work well in predicting the binding affinities of ligands because the binding sites of proteins are not flexible. If the environment for a property is rigid (e.g., the binding sites of proteins), the model for predicting the property should consider H-bond interactions individually, rather than the overall H-bond interactions. In our further study, we will explore how to develop models for the properties related to rigid environments, which can be used to develop scoring functions for predicting protein-ligand binding affinities and develop QSAR models for screening databases of ligands. Furthermore, in this study, we demonstrated to effectiveness of the LFER for predicting the properties of neutral organic compounds. If a dataset contains ionizable compounds, it will be necessary to include molecular descriptors for the ionized forms. Although several approaches currently exist for considering the effects of ionization on various molecular properties (Li et al., 2006; Zhang et al., 2017), our future work with involve attempts to adapt the proposed LFER for use in these situations.