Chem. Res. Toxicol. 1988,1, 123-127
123
A Structure-Activity Relationship Study of Organophosphorus Compounds R. H. Rohrbaugh and P. C. Jurs* 152 Davey Laboratory, Chemistry Department, The Pennsylvania State University, University Park, Pennsylvania 16802
W. P. Ashman, E. G. Davis, and J. H. Lewis Research Directorate, Chemical Research Development and Engineering Center, Aberdeen Proving Ground, Maryland 21010 Received January 11, 1988
Organophosphorus compounds have been shown to exhibit toxic behavior as insecticides, pesticides, and mammalicides. Soman and 21 related compounds were studied for possible structure-activity relationships. Computer-aided methods were used to generate a linear expression relating the activity (In [l/LDm],rabbit I.V.) of the compounds to three structure-based descriptors (R = 0.96). Principal components regression and jackknife analysis were performed to assess the stability of the model.
Background Organophosphorus compounds have been used extensively as insecticides, pesticides, and mammalicides. These compounds act as inhibitors in the hydrolysis of acetylcholine (ACh). Acetylcholine is an important neurotransmitter which is commonly found in the neuromuscular junctions. The presence of ACh is extremely important to the functioning of living organisms, especially in respiration. The organophosphorus compounds can act by covalently binding with the acetylcholinesterase (AchE) enzyme and thus preventing the hydrolysis of acetylcholine to acetate and choline (Figure 1). This disruption of the acetylcholinergic system can lead to partial or complete paralysis, as well as death by asphyxiation. Quantitative structure-activity relationship studies (QSAR) have been used to describe the effects of organophosphorus compounds on the cholinergic system. Magar used QSAR methods to identify parameters which model the inhibition of organophosphorus compounds on myoneural junctions (1)and organophosphorus compound toxicity (2). Simon and co-workers have used QSAR methods to relate toxicity of organophosphoruscompounds with physicochemical parameters (3), and Johnson reported on the use of quantum chemical parameters to predict the relationship of organophosphorus compound structure to acetylcholinesterase inhibition (4). Effective medical treatment for organophosphorus poisoning usually consists of an oxime to reactivate the inhibited acetylcholinesterase enzyme. Soman (I) is an O2 1
I13
H
I I
6
CH3-P-O-C-C-CH3
I
F4
CH3
CH3
/
/
CH3
I
extremely toxic organophosphorus compound that is very resistant to current organophosphorus therapeutic treatment (5). In an effort to determine molecular features of importance for designing new and effective therapies against soman poisoning, a QSAR study was performed 0893-228~/88/2701-0123$01.50/0
to identify parameters related to the LDbotoxicity in a series of soman type compounds. A number of physicochemical and steric parameters were calculated and evaluated. A QSAR model relating molecular shape and molecular connectivity has been developed that correlates highly with toxicity.
Materials and Methods Data Set. The data set consisted of toxicity data for the 22 compounds in Table I. The measured toxicity consisted of rabbit intravenous LDM values expressed as In (l/pmol/kg) (6). The structures and toxicity data were entered into computer files for analysis using the ADAPT software system (7). Molecular Modeling. The molecular modeling analysis and display system (8) was used to perform Allinger's MM2 (9, 10) calculations to determine the minimum energy conformation of each molecule. The resultant three-dimensional optimized conformations (11)were used for the generation of three dimensional molecular descriptors. The parameter set required for structure geometry optimization supplied with the MM2 program was used. If the required parameters were not available within the MM2 program, ab initio calculations using Gaussian82 (12) were used to generate the needed parameters and/or parameters were provided by Leonard (13) and Hopfinger (14). Physicochemical Descriptors. The lipophilic constant r (fragmental octanol/water partition coefficient) was calculated by using the MedChem (15) Project's CLOGP program. The r values for the phosphorus functional groups in the molecule were provided by Famini (16). The lipophilic constant is often useful for approximating the partitioning abilities of a solute in biological systems. The molar refraction values were calculated by using the program CMR that is contained in the MedChem Project. Molar refraction is related to the polarizability of the molecule and is considered a bulk property descriptor. Topological Descriptors. Several topological descriptors were generated for use in this study. These included molecular connectivities (x),kappa indices ( K ) , and a branching indicator. Molecular connectivity indices were originally developed by RandiE (17) and later modified by Kier and Hall (18). These indices are based on a graph theoretical treatment of the molecular topology of the compounds, and encode information about the branching and size of the molecules. The general equation for 1988 American Chemical Society
124 Chem. Res. Toxicol., Vol. 1, No. 2, 1988
fI
Enzyme -Ser - O H
Rohrbaugh et al. Table I. Structures and Activities
+ H3C-C-O-CH2-
0
CH,- N-(CH,),
II I
CH3-P-OX F
+ H O - C H r CH2- N7CH313 0
x 1
II Enzyme -Sei -0-C-CHj
i H2O
2
C
0
Enzyme - S e i -OH
-cccI
2.830
-ccc
2.645
1
3
4
+ CH,-C-O
I C
C
2.639
1
-ccc I
Figure 1. Hydrolysis of acetylcholine (ACh) by acetylcholinesterase (AChE).
5
calculating molecular connectivities of the nth order is given in eq 1,where "xt is the nth order term of type t (t = path, cluster,
6
7 s=l
2.919
-ctc II cc
H+
II
activity
-ccc I c-c
2.564
F
2.386
-ccc 1 ccc -cccc I
2.354
i
path-cluster, or chain), n, is the number of connected subgraphs of type t, m is the number of edges, and 6 is the vertex valence. The following valence molecular connectivities were calculated in the present study: O x , ' x , 2x, 3xp,3xpe,4xp,and 4xqc. Six kappa indices were also calculated. These indices were developed by Kier (19,20) and encode topological shape using a graph theoretical approach. Each structure is depicted as a graph consisting of nodes (atoms) and edges (bonds). For the fmt three kappa indices, no distinction is made between different atom and bond types. The 1~ index is based on the number of paths of length one (or the number of one bond fragments) found in the molecule. Similarly, the % and ' K indices are based on the number of paths of length two and three. For each index, the number of paths of the appropriate length is taken in relation to the minimum and maximum number of paths of the same length possible for a graph containing the same number of nodes. In general, the nth order kappa index, "K, for a molecule with N atoms is given by eq 2, where C is a constant, "P,, is the max-
8 9
10
2.313
-cc
2.235
2.302 2.273
I
11 12 13
14 15 16
imum number of paths of length n for N nodes, "Pmhis the minimum number of paths of length n for N nodes, and "Pi is the actual number of paths of length n for the structure. Three were also used which were designed modified indices (1K,,2K,,3K,) by Kier to take atom and bond types into consideration (21). An indicator variable which encoded the presence of branching was also used. If branching from the main chain occurs in the substituent X (Table I) of each structure, a value of 1 is assigned to the variable, otherwise, a value of 0 is assigned. Shape Descriptors. Molecular Shape Analysis (MSA). Molecular shape analysis (MSA) was performed to encode the shape related information in the molecular structure of the compounds. MSA is achieved by superimposing each structure in the data set with a selected reference structure. Typically, the reference structure chosen is that compound which exhibits the highest observed activity. For the study at hand, soman (I) was chosen since it has the highest toxicity value (2.919). Hopfinger's MSA program (22) was used to calculate the volume of each compound, the MSA overlap volume (MSA OV), and the MSA nonoverlap volume (MSA NOV). The MSA overlap volume is calculated by taking the soman molecule (I) and overlapping atoms 1-36 with the corresponding atoms 1-3-5 of each of the other compounds studied. The MSA program computes the resultant overlap volume. The MSA nonoverlap volume is calculated by subtracting the overlap volume of each individual compound from its total volume. Another descriptor evaluated, the compound's overlap volume difference, is determined by subtracting the compound's MSA overlap volume from soman's volume.
-cccc I 1 c c c-c, -c: c c-c' F-c, -c, c c-c'
17 18 19
20 21 22
I
C
-cccI c-c, -c \ Ic-c c-c -ccc C
I
-cccc I 1 c c -cccc -cc -C C I
-cccc I ccc -cccccc -cccccc I cc -ccccccc I C
2.096 1.938 1.726 1.590 1.127
1.030 0.983 0.978 0.798 0.629 0.511
Shadow Areas. Shadow areas (23,24)were also calculated for the 22 compounds in an effort to encode shape information. Shadow areas consist of the areas of two-dimensional projections of a molecule onto the planes defined by the X , Y , and 2 axes. In order to make comparisons, the compounds must all be aligned in some unambiguous manner prior to calculation of the shadow areas. AlJ compounds studied were oriented such that atoms 1, 2, and 3, were placed in the X-Y plane. For each compound, 12 parameters related to these areas were calculated. Shadow areas 1-3 (Sl-S3) consist of the areas of the projections onto the X-Y, X-2, and Y-2 planes, respectively (Figure 2). Shadow areas 4-6 (S4-S6) are normalized for the size of the molecule. For this study six new parameters were developed. The new shadow areas are utilized in a comparative manner, similar to MSA. For each compound, the three projections were overlapped with the three corresponding projections of a standard compound, and the overlapped (S7-S9) and nonoverlapped (SlO-Sl2) areas
Chem. Res. Toxicol., Vol. 1, No. 2, 1988 125
Structure-Activity Study of Organophosphorus Compounds
Table 111. Multiple Linear Regression Resultsa 95% CI VIF descriptor coeff.
3.580
0
o
> H
=
?I
s
0 ~
W
1.958
-
r
0
r
n
o
ff
3 -.
Figure 2. Shadow areas S1, S2, and S3 for soman (I).
a a
0
0
mu
3
0.400
I ,950
3.508
MEASURED A C T I V I T Y
Figure 4. Plot of estimated versus measured activity wing normal regression model. several Variables. Several linear regression techniques were applied to the data in order to determine the best subset of molecular descriptors which could be used to describe the activity of the organophosphorous compounds. These techniques included stepwise multiple linear regression (25,261 and leaps and bounds (27).
Figure 3. Overlap (crosshatch) and nonoverlap (/// hatch) shadow areas (S7 and SlO) for compound 20 (top). Table 11. Correlation Matrix for Final Descriptor Pool MSA OV 1.000 0.755 1.000 ?%/MI3 -0.231 -0.280 1.000 0.263 0.468 LOO0 S8 0.028 0.497 0.392 0.475 1.000 s11 0.615 for each projection were calculated. As with the MSA analysis, soman (I) was used as the standard structure. The overlap and nonoverlap areas are illustrated for compound 20 in Figure 3. Moments of Inertia. The three principal momenta of inertia were calculated for each compound using the moment of inertia tensor method. For each compound, the largest (MIl), second largest (MI2), and smallest (MI3) moment of inertia were calculated, as well as their ratios. Initial Descriptor Selection. Thirty- eight descriptors, as defined previously, were calculated for consideration. Initial descriptor selection involved deletion of those variables which were involved in high simple and multiple correlations (R > 0.9) with each other or in combination. After the initial deletion step, only five variables remained for analysis (Table 11). The highest remaining simple correlation was R = 0.755 between MSA overlap volume and 4xpcmolecular connectivity. The highest remaining multiple correlation was R = 0.864 between MSA overlap volume and 4xpc,S11, and S8. Regression Analysis. Multiple linear regression (MLR) can be used to investigate the presence of linear relationships between
Results and Discussion The regression analysis results are given in Table 111. The equation shows a strong linear relationship between the activity values and the MSA overlap volume, ,511,and 4xpc. The multiple R of the equation is 0.96, and the standard error is 0.24 (approximately 10% of the mean activity). Figure 4 shows a plot of the estimated activities versus the experimentally measured activities. As previously stated, collinearity was still present to a certain extent in the data when regression analysis was performed. In order to evaluate the effect of removing one of the two most highly correlated descriptors (MSA overlap volume and 4xpc),two additional regressions were performed. In each case, the resultant equation was measurably inferior (based on multiple R and standard errors) to the original equation. As a result, it was decided to retain both descriptors and to try to address the collinearity problem in another manner. In an effort to determine the effects of the collinearity on the quality of the estimated coefficients, the variance inflation factors (VIF) were calculated for each term in the original model by using MINITAB. These values are reported in the last column in Table 111. The values indicate that the variances associated with the estimation of the coefficients are 180%, 60%, and 130%, respectively, for the three terms. On the basis of this information, quantitative interpretation of the meaning of the variables selected is not justified. However, the presence of the higher VIFs (as compared to completely orthogonal variables) does not necessarily diminish the predictive abilities of the model (28). In order to verify this, an alternative regression study was performed by using princippl components.
126 Chem. Res. Toxicol., Vol. 1, No. 2, 1988
PC 1 2
3 4 5
eigenvalue 400.81 39.26 7.80 4.62 0.06
Rohrbaugh et al.
Table IV. Principal Components Analysis of Final Descriptor Pool 1oadings cum. var. MSA OV 'XPC MIl/MB S8 -0.0174 0.0342 -0.0113 88.56 -0.976 -0.00219 97.24 0.157 -0.442 -0.456 98.96 -0.0644 -0.0502 0.244 -0.888 -0.0385 99.98 0.138 0.861 0.0205 100.00 0.0146 -0.998 -0.0451 0.0450
Table V. PrinciDal ComDonents Regression Resultsa descriptor coeff. 95% CI PC1 -0.025 f0.0026 PC2 0.080 f0.0082 PC5 -0.70 f0.20 PC4 0.056 f0.024 intercept -1.3 O n
~~
s11 -0.216 -0.756 0.380 -0.487 0.00505
Table VI. Jackknife Analysis Results normal DrinciDal comDonents R" 0.93 0.92 1.82 X 10" mean residual -3.27 x 10-3 R* 0.96 (0.32%) 0.96 (0.38%) SC 0.236 (3.9%) 0.235 (5.0%) % Deviation in Regression Coefficients
= 22; R = 0.96; s = 0.24.
principal components
normal
3.580
intercept term 1 term 2 term 3 term 4
MSA OV
s11
4xpe
10.9% 2.5% 3.5% 9.4%
PC1 PC2 PC5 PC4
11.8% 5.5% 3.8% 13.5% 14.9%
P
Correlation between activities and jackknifed y-hats. bAverage R for all jackknifed regressions (% SD in parentheses). CAverages for all jackknifed regressions (% SD in parentheses).
H H !4
i3
I ,950
U H
B
0,480 0,400
1.958
3.500
MEASURED A C T I V I T Y
Figure 5. Plot of estimated versus measured activity using principal Components regression model.
Principal Components Regression. Principal components analysis (PCA) (29) was performed on the five variables in the final descriptor pool. The eigenvalues and cumulative variances for the five principal components (PCs) are given in Table IV. Also in Table IV are the individual loadings of each descriptor into each PC. By definition, the five principal components contain all of the information present in the original variables and yet are completely orthogonal. As a result, there should be no adverse effects on the estimation of the coefficients of the variables due to collinearity. Performing the regression analysis on the five principal components yielded a fourvariable equation (Table V) with a multiple R of 0.96 and a standard error of 0.24. These values are identical with those obtained with the original variables. Figure 5 gives a plot of the predicted versus experimentally measured activities using the principal Components regression results. Comparison of Equations. The four PCs selected on the regression analysis (in order of importance) were PC1, PC2, PC5, and PC4. By examining the loadings in Table IV, it is evident that PC1 is predominantly loaded with MSA overlap volume (-0.976),PC2 with S11(-0.756), PC5 with 4xpc(-0.998), and PC4 with MIl/MI3 (0.861). The first three PCs selected are heavily loaded with the same descriptors which were selected in the normal regression analysis, as well as being in the same relative order of importance. This observation, coupled with the identical statistics (R and s) of the two equations, suggests that both the normal regression and the PC regression are equally
effective in describing the activity of the compounds. The next step in comparing the two models was to examine both their validity, as well as the stability of the models' coefficients. This was achieved by using jackknife analysis (25). In jackknife analysis, each of the n compounds is held out of the analysis in turn. A regression model is built on the basis of the remaining n - 1 compounds. The activity of the remaining compound is then estimated on the basis of the derived model. In addition to giving a rough idea of the predictive ability of a model, jackknife analysis can also be used to evaluate the stability of the models' coefficients with respect to the individual observations. Table VI gives a comparison of the results of the jackknife analysis for the two regressions. The correlation between the jackknifed estimated activities and the measured activities was 0.93 for the normal regression and 0.92 for the PC regression. Comparing the mean jackknifed residuals for the two regressions shows a value of -3.27 X for the normal regression and 1.82 X for the PC regression. A comparison of the averages of the individual jackknifed regression statistics shows that the two models are essentially equivalent; although the normal regression has slightly better statistics than the PC regression. Finally, through examination of the percent deviation in the coefficients of each term during the jackknife analysis, it is evident that the PC regression coefficients vary to a higher degree than do the normal regression coefficients. This would indicate that the normal regression is slightly more stable with respect to the individual observations than is the PC regression. Ideally, prediction studies should be performed to assess the usefulness of each model. However, because of the limited size of the data set, no compounds were set aside for use in predictive studies. Despite this fact, the equations generated may still be used to estimate the activities of untested compounds. However, the proposal of new active compounds was not within the scope of the present work. Interpretation of Molecular Descriptors. Upon examination of the descriptors in the linear models, it
Structure-Actiuity Study of Organophosphorus Compounds
becomes evident that the conformation or three-dimensional geometry of the compounds is important in describing the observed activities. The MSA overlap volume and the shadow area parameter (Sll) are both shape related parameters. The 4xpcterm is, to a certain degree, a measure of the presence of midchain branching. For this series of compounds, branching tends to make the compounds more spherical in shape, whereas the absence of branching will tend to make the compounds more linear. In this respect, the 4xpcterm is also an indirect measure of shape. The importance of the two comparative terms (MSA overlap volume and S11)suggests the existence of an optimal shape/conformation, which most likely mimics the receptor site on the AChE enzyme. The active site of the AChE has been determined to consist of an anionic and an esteratic region (30-32). In order for a potential inhibitor to exhibit activity, the shape / conformation must be such that the necessary interactions (Coulombic, hydrophobic, etc.) can occur.
Conclusions Two linear models were generated relating the LDW activity of 22 organophosphorous compounds to several structure-based descriptors. Both equations were statistically significant with multiple correlation coefficient values of 0.96 and standard errors of 0.24. The equation generated by using untransformed variables contained a fair degree of collinearity and yet was more stable than the equation generated by using principal components. This suggests that although the collinearity present in the untransformed variables is such that quantitative interpretation is not feasible, the predictive ability and stability of the model does not suffer. Finally, the descriptors chosen for the linear model suggests the importance of the three-dimensional conformation of the compounds in determining the activities of the organophosphorus compounds. Acknowledgment. This work was supported by the National Science Foundation under Grant CHE-8202620. The PRIME 750 computer was purchased with partial financial support of the National Science Foundation. R.H.R. also thanks the American Chemical Society Analytical Division Fellowship Program for financial support. Registry No. 1,96-64-0; 2,6154-51-4; 3,352-52-3; 4,372-62-3; 5,66348-71-8; 6,113548-85-9; 7,761-93-3; 8,352-53-4; 9,329-99-7; 10, 113548-86-0; 11, 107-44-8; 12, 2053-81-8; 13,113548-87-1; 14, 763-14-4; 15,30593-65-8; 16, 352-63-6; 17,673-97-2; 18,353-88-8; 19, 113548-88-2; 20, 113548-89-3; 21, 458-71-9; 22, 22925-97-9.
References (1) Mager, P., and Seese, A. (1982) “Organophosphorus pesticides and mammalicides”. Zool. Jahrb., Abt. Anat. Ontog. Tiere 107, 46-70. (2) Mager, P. (1981) “Grouped observations in quantitative structure-toxicity relationships of organophosphorus poisons”. Pharmazie H 36, 513-514. (3) Simon, Z. (1980) “Steric fit in quantitative structure-activity relations”. Lecture Notes in Chemistry (Balaban, A,, Chiriac, A., Motoc, I., and Simon, Z., Eds.) Vol. 15, pp 67-73, Springer-Verlag, Berlin. (4) Johnson, H. Kenley, R. A., Rynard, C., and Golub, M. A. (1985) -QSAFt for cholinesterase inhibition by organophosphorus esters and CND0/2 calculations for organophosphorusester hydrolysis”.
Chem. Res. Toxicol., Vol. 1, No. 2, 1988 127 Quant. Struct. Act. Relat. Pharmacol., Chem., Biol. 4, 172-180. (5) Boskovic, B. (1981) “The treatment of soman poisoning and ita perspectives”. Fundam. Appl. Toxicol. 1, 203-213. (6) Data provided by Toxicology Division, Research Directorate, Chemical Research Development and Engineering Center, Aberdeen Proving Ground, MD. (7) Stuper, A. J., Brugger, W. E., and Jurs, P. C. (1979) Computer Assisted Studies of Chemical Structure and Biological Function, Wiley-Interscience, New York. (8) Leonard, J. J. (1986) CRDEC-TR-86039. A User’s Guide to the Molecular Modeling, Analysis and Display System (MMADS), Chemical Research Development and Engineering Center, Aberdeen Proving Ground, MD. (9) Allinger, N. L., and Yuh,Y. H. (1980) QCPE 13,395. (10) Allinger, N. L. (1973) “Conformational analysis. 130. MM2. A hydrocarbon force field utilizing V1 and V2 torsional terms”. J. Am. Chem. SOC. 95, 8005-8025. (11) Ashman, W. P., Research Directorate, Chemical Research Development and Engineering Center, Aberdeen Proving Ground, MD (unpublished data). (12) Binckley, J. S., Whiteside, R. A, Raghavachari, K., Seeger, R., Defrees, D. J., Schlegel, H. B., Frisch, M. J., Pople, J. A., and Kahn, L. R. (1982) Gaussian82 Release A, Carnegie Mellon University, Pittsburgh, PA. (13) Leonard, J. M., Research Directorate, Chemical Research Development and Engineering Center, Aberdeen Proving Ground, MD (unpublished data). (14) Hopfinger, A. J. (unpublished data). (15) MedChem Software, Release 3.33. Medicinal Chemistry Project, Pomona College, Claremont, CA (1985). (16) Famini, G., Research Directorate, Chemical Research Development and Engineering Center, Aberdeen Proving Ground, MD (unpublished data). (17) RandiC, M. (1975) ”On characterization of molecular 97,6609-6615. branching”. J. Am. Chem. SOC. (18) Kier, L. B., and Hall, L. H. (1986) Molecular Connectiuity in Structure-Activity Analysis”. Wiley, New York. (19) Kier, L. B. (1985) “A shape index from molecular graphs”. Quant. Struct.-Act. Relat. Pharmacol., Chem., Biol. 4, 109. (20) Kier, L. B. (1986) “Shape indexes of orders one and three from molecular graphs”. Quant. Struct.-Act. Relat. Pharmacol., Chem., Biol. 5, 1-7. (21) Kier, L. B. (1986) “Distinguishing atom differences in a molecular graph shape index”. Quant. Struct.-Act. Relat. Pharmacol., Chem., Biol. 5, 7-12. (22) Hopfinger, A. J. (1980) “A QSAR investigation of dihydrofolate reductase inhibition by Baker triazines based upon molecular shape analysis”. J. Am. Chem. SOC.102, 7196-7206. (23) Rohrbaugh, R. H., and Jurs, P. C. (1987) “Molecularshape and the prediction of high-performance liquid chromatographic retention indexes of polycyclic aromatic hydrocarbons”. Anal. Chem. 59, 1048-1054. (24) Rohrbaugh, R. H., and Jurs, P. C. (1987) “Descriptions of molecular shape applied to structure-activity and structureproperty relationship studies”. Anal. Chim. Acta 199, 99-109. (25) Draper, N. R., and Smith, H. (1981) Applied Regression Analysis, 2nd ed., Wiley, New York. (25) Belsley, D. A., Kuh, E., and Welsh, R. E. (1980) Regression Diagnostics: Identifying Influential Data and Sources of Collinearity, Wiley, New York. (27) Furnival, G. M., and Wilson, R. W. (1974) “Regression by leaps and bounds”. Technometrics 16, 499. (28) Neter, J., Wasserman, W., and Kutner, M. H. (1985) Applied Linear Statistical Models, 2nd ed., pp 335-393, Richard D. Irwin Inc., Homewood, IL. (29) Johnson, R. A., and Wichern, D. W. (1982) Applied Multiuariate Statistical Analysis, Chapter 8, Prentice-Hall, Englewood Cliffs, NJ. (30) Heath, D. F. (1961) Organophosphorus Poisons. Anticholinesterase and Related Compounds, Pergamon, New York. (31) O’Brien, R. D. (1964) Toxic Phosphorus Esters, Academic, New York. (32) Englehard, N., Prchal, K., and Nenner, M. (1967) “Acetylcholinesterase”. Angeur. Chem., Znt. Ed. Engl. 6, 615.
