The forensic medical examination of poisoning cases with psychoactive substances is carried out using targeted analytical methods, such as gas chromatography–mass spectrometry (GC-MS) and liquid chromatography-tandem mass spectrometry (HPLC-MS/MS), to detect known xenobiotics and their metabolites [1]. Untargeted testing using high-resolution liquid chromatography-mass spectrometry (HPLC-HRMS) technologies can be used to identify the chemical structure of new substances appearing in illicit traffic. Substances identified in biological samples can pose a significant threat [2–5]. Although the safety of a particular drug can be measured by determining its therapeutic index, calculated as the ratio of the drug maximum dose that does not exhibit toxicity to the dose that provides the desired effect (LD50/ED50), this approach cannot be used for substances detected during forensic medical examination. In addition, information about the hazardous pharmacological activity of xenobiotics may frequently be either fragmentary and contradictory or absent altogether, for a number of reasons:

The above reasons hamper the medicobiological assessment of new psychoactive substances, although there has been a growth in the number of cases of illicit trafficking both in Russia and globally [4][8–10]. In order to be able to predict the potential hazard of CNS-active xenobiotics and other chemical substances, criteria for a straightforward determination of their toxic effect are required.

It is common knowledge that the toxicity of many narcotic analgesics is associated with opioid-induced respiratory depression, which result in the patient’s death in the absence of proper treatment [11]. The central generator of breathing patterns in the brainstem is the preBötzinger Complex and the Kolliker–Fuse nucleus. In these areas of the brain, only 70–140 neurons are involved in responding to an increase in carbon dioxide (pCO2) levels and a decrease in blood oxygenation (pO2), which is necessary for the reflex mechanism of respiration [12]. Accordingly, the mechanism of regulation of respiratory activity is extremely vulnerable to substances capable of selectively affecting these targets in the brain. It can be assumed that the health hazard of CNS-active xenobiotics, expressed in terms of strength and speed, is mediated, on the one hand, by factors determining the degree of their effect on target receptors in the neurons of these brain regions, and, on the other, by factors determining the transfer of xenobiotics from the central bloodstream to their biological targets in the brain.

The strength of narcotic analgesics is often associated with the inhibition constants of opioid receptors [13]. However, under standardized conditions, morphine and fentanyl, which differ significantly in the strength of their analgesic effects, possess comparable inhibition constants (Ki) of the μ-opioid receptors (MOR) (Table 1) [5][7].

Upon closer examination, no direct relationship exists between the analgesic activity of opioids and the strength of their binding to the receptor, which was earlier confirmed by numerous examples [14]. Thus, according to Boström et al., oxycodone, which binds to opioid receptors more than 20 times less strongly than morphine (Table 1), exhibits a 1.8-fold greater analgesic activity. This effect can be explained by its higher concentration in the target brain tissues than in the blood. Morphine, on the contrary, has an achievable concentration in the brain of rats several times lower than in the blood [15]. Thus, the binding of small molecules to opioid receptors in most cases is the determining condition for manifestation of analgesic activity, although being not the main factor in the strength and speed of its manifestation (similar to the classical works of N.V. Lazarev on the manifestation of narcotic properties in hydrocarbons (non-electrolytes)) [16]. The hazardous effects of opioid analgesics are also related to their properties, which determine the characteristics of entry into target tissues, organs, and targets in the central nervous system (CNS) through the blood-brain barrier (BBB).

The BBB acts as a filter through which nutrients flow from the bloodstream to the brain and in the opposite direction, while the waste products of the nervous tissue are removed. The BBB protects the brain from microorganisms, toxins, cellular and humoral factors of the immune system, and xenobiotics circulating in the blood. Drugs acting on CNS targets must have the ability to penetrate such a biological barrier. Currently, three main ways of transporting small molecules to brain tissues are being considered, including unidirectional penetration into the brain through passive and facilitated diffusion, due to a concentration gradient of substances and requiring no additional energy; active transport, requiring energy (ATP) for the transport of molecules against a concentration gradient; passive diffusion for moderately lipophilic medications and active transport mechanisms for penetration through the BBB of predominantly polar molecules [17][18].

The membranes separating the cells of the BBB layers act as a channel for the diffusion of organic molecules through the BBB cellular layers. In order to pass through such a channel, a molecule dissolved in the blood must be transformed during a multi-stage physicochemical process of desolvation and charge loss for a non-covalent interaction with the phospholipid bilayer. This process is similar to that of ligand-receptor interaction, although with lower steric requirements. Therefore, for penetration into the brain, small molecules must exhibit a set of specific properties, such as optimal geometric size, ionization properties, flexibility, etc. Thus, lipophilicity is a parameter that correlates well with the analgesic power of anesthetics and analgesics. This parameter, expressed in the form of LogP [16], is now considered as a composite descriptor describing the contribution of steric intermolecular interactions and that of formed hydrogen bonds [17]. Descriptors of hydrogen bond formation are properties of molecules, such as polarizability, polar surface area, the number of donors and acceptors of hydrogen bonds, or heteroatoms capable of their formation. Polar compounds with a high potential for hydrogen bonding, e.g., peptides, cannot easily penetrate the BBB. Polarizability, underlying a number of physical properties of substances, including surface tension and solubility, is characterized by a dipole moment. The polar surface area is defined as the sum of the surfaces of all polar atoms or molecules, primarily oxygen and nitrogen. Medications acting on the CNS possess smaller polar surface areas (PSA) than other classes of drugs, usually no more than 90 A2. Molecules with a polar surface area of more than 140 A2 do not readily penetrate cell membranes.

Geometric and steric factors expressed by the molecular weight, molecular volume, and flexibility of the molecule, are important for diffusion through biological membranes. For effective penetration through the BBB, the molecular weight should not exceed 400 Da. The molecular volume as a function of the molecular weight and structure takes all the conformations available to the molecule into account. The presence of more than ten rotating connections in structures correlates with a decrease in their bioavailability and CNS activity. The limited flexibility and compactness of molecules with fewer polar groups on their surface capable of functioning as donors and acceptors of hydrogen bonds is an advantage for substances acting on the CNS; therefore, the parameters describing these properties are significant factors related to the CNS activity of xenobiotics.

Most CNS-active substances have a nitrogen atom capable of ionization in the blood plasma. The resulting charge of the protonated molecule negatively affects its ability to diffuse through the dielectric layers of phospholipid membranes, compared to neutral molecules [18]. However. the presence of a positive charge on the nitrogen atom is considered an important factor for the electrostatic interaction of substances with the key amino acids of target receptors [19]. The molecule ionization degree in the buffer system of the blood and cerebrospinal fluid is a dualistic factor related both to the rate of entry of active molecules into CNS tissues and influencing their ability to bind to CNS targets [20, 21]. Thus, the severity of the toxic effect of CNS-active substances, due to the effectiveness of their entry into the brain, is determined by the sum of physicochemical factors related to their chemical structure. These factors can be experimentally evaluated or calculated.

This work is aimed at studying the relationship between the physicochemical properties of one class of CNS-active substances and their health hazard potential using mathematical analysis. On this basis, criteria for preliminary hazard assessment of narcotic drugs are substantiated.

Table 2 presents a sample of model objects used in the study selected from the known structures of narcotic analgesics. The objects under consideration were divided into three groups:

The physicochemical properties that can be measured without laborious experiments and reference materials of controlled substances were considered as potential descriptors of the hazard rate of CNS-active xenobiotics (see Table 3). The molecular weight of substances was determined by GC-MS; the distribution constants “octanol/water” and acid-base dissociation correlated with retention times were established by HPLC [21][23]. In cases where information about the structural formula of a xenobiotic can be found in mass spectrometry databases, computer simulation methods using accessible and reliable calculation algorithms are used to determine its physicochemical properties [24].

The values of the selected physicochemical parameters for the training sample of substances were calculated using the ACD/Percepta software [24]. The methods of descriptive statistics and the method of linear discriminant analysis using the Statistica 6.0 statistical analysis application were used for mathematical data analysis. We previously showed the effectiveness of multidimensional statistics (chemometry) methods, in particular linear discriminant analysis and multiple regression for calculating toxic properties of refrigerants based on the physicochemical parameters of molecules [25], for assessing the degree of exposure to xenobiotics and noninvasive diagnosis of a number of diseases based on the spectral characteristics of biological samples1 [26].

The generated sample contained 53 anonymized records of the calculated values of the physicochemical properties of the model CNS-active substances described in the literature (Table 2). The structures of the analyzed model preparations were used for calculating the selected physicochemical quantities, the results of which are shown in Table 4.

Table 5 presents the average values of the physicochemical descriptors used to construct regression models and to determine the limits of their applicability.

The average values of the physicochemical quantities presented in Table 5 show that the narcotic substances in Group 2 demonstrate the highest lipophilicity (LogP). At the same time, the ratio of the average lipophilicity values of the drug groups in water and blood plasma — LogP and LogD7.4 — shows that the lipophilicity of Groups 1 and 2 significantly decreases at physiological pH. Conversely, the highly active substances belonging to Group 3 and having a lower value of the basic ionization constant demonstrate high lipophilicity not only in water, but also in blood plasma. This emphasizes the relationship between lipophilicity and the ability of substances to enter CNS tissues during redistribution in the blood-brain system.

The physicochemical criteria that can be used for hazard assessment of CNS-active substances were established based on paired correlation coefficients between the physicochemical descriptors shown in Table 6. It can be seen that most descriptors are weakly correlated with one another. The strongest correlation is observed between the LogD7.4 and LogP values since these are highly similar parameters in nature. The difference between them is related to the pH values at which they are determined, being due to the degree of ionization of the molecules of substances at these pH values. This is indicated by a rather strong correlation (R = 0.97) between the calculated value (LogD7.4 – LogP) and the pKb value.

A strong correlation is also observed between the values of MW and Polar. This is likely to be due to a proportional increase in the integral polarizability of the molecule, associated with an increase in the number of polarizable fragments with increasing mass. An analysis of the obtained correlations between the physicochemical properties in the training groups of CNS-active substances showed that the values considered, either by themselves or in pairs, do not allow for an unambiguous distribution of substances to one of the three groups of hazardous substances.

A discriminant analysis was used to analyze the entire data set of physicochemical properties of model substances in order to identify their relationship with the degree of their health hazard. Discriminant analysis, as a branch of multidimensional statistical analysis, includes statistical methods for classifying multidimensional observations in a situation where the researcher has a priori so-called training samples (classification with learning). Discriminant analysis makes it possible to classify an object based on the measurement of various characteristics (features, descriptors), i.e., permitting their assignment to one of several groups (classes). In this study, discriminant analysis was used to determine the differences of aggregates in the mean of a variable (or a linear combination of variables) for subsequent use of this variable or a linear combination of variables as a criterion for the membership of new members to a particular group.

The resulting data set was subjected to a linear discriminant analysis (LDA) procedure with step-by-step inclusion of variables. The grouping variable was the “Group” category, with all other descriptors being used as independent variables. Initially, the data array was divided into two parts: one part was used as a training part to derive a mathematical classification model, while the second part acted as a control part to verify the resulting model. The last three compounds in the list of each group were used as controls.

In the process of LDA with step-by-step inclusion of variables, discriminant functions were determined. At each step, we analyzed all the variables to find the one that made the greatest contribution to the difference between the aggregates. This variable was included in the model at this step, and then the transition to the next step was carried out. The step-by-step variable inclusion mode allowed the minimum of variables to be involved in the analysis. The discriminant functions, expressing the maximum heterogeneity of groups among themselves, were linear combinations of variables optimized such that to discriminate groups among themselves most effectively. Since all discriminant functions were orthogonal, they collectively formed the phase space of discriminant functions, i.e., an n-dimensional Cartesian coordinate system, where n is the number of statistically significant discriminant functions.

The discriminant functions were uniform for the data set involved in LDA; however, the values of the roots of each of the discriminant functions were strictly specific for each substance. The roots of discriminant functions for each substance were calculated by matrix multiplication of a vector of descriptors of this substance by a matrix of coefficients of discriminant functions. Thus, the original data matrix was transformed into a matrix containing values of the roots of discriminant functions instead of descriptors. The values of the roots, unique to each junction, were essentially the coordinates of the points of the corresponding junctions in the phase space of the discriminant functions. The average values of these roots for each group determined the coordinates of the group centers (the so-called group centroids). It was the values of the roots of discriminant functions or the coordinates of objects (substances) in the phase space of discriminant functions that made it possible to classify objects, i.e., to correlate a specific substance with any of the groups by the proximity of a point to the centroid of the group.

The conducted LDA of the training sample of the data matrix with descriptors for three groups of substances produced two statistically significant discriminant functions dF1 and dF2. These functions were fully descriptive the data array, of which the dF1 function describes 88% of the information contained in the data presented in Table 7.

Thus, the obtained discriminant functions are the following linear combinations of descriptor values:

dF1 = 0,051×MW + 2.599×LogP – 1.943×LogD7.4 – 0.643×Polar – 0.032×TPSA + 6.638

dF2 = 0.028×MW – 1.376×LogP – 0.845×LogD7.4 – 0.359×Polar + 7.010

The standardized coefficients for each descriptor for the dF1 function, which are relatively close to each other in absolute value and far from zero, indicate the important contribution of all the considered physicochemical indicators to the discrimination performed by this most statistically significant discriminant function. In principle, the same applies to the dF2 function, with the exception for the TPSA (topological area of the polar surface) indicator, which proved to be negligible.

For all model substances from Table 4, including the control ones, the values of the roots of the discriminant functions dF1 and dF2 were calculated. The root values, unique for each junction, are the coordinates of the points of the corresponding junctions in the phase space of discriminant functions, which is a flat coordinate system with the dF1 and dF2 axes. Figure 1 demonstrates the location of points corresponding to specific samples of substances from the three groups in the phase space of the dF1 and dF2 functions.

It can be seen from Fig. 1 that all groups are well separated from each other; all control points were determined in their groups.

For all three groups of model substances, the average values of the roots of the dF1 and dF2 discriminant functions were found, which are the coordinates of the centroids of the groups in the phase space of the dF1 and dF2 functions (Tables 8 and 9). The obtained discriminant functions, according to the data of the classification matrix (Table 9), possess a high discriminating power, allowing a fairly clear division of the groups of substances under consideration among themselves. Table 9 shows the LDA classification results for both the substances in the training sample and for all substances, including the control ones.

In general, the substances were classified correctly, with the exception for substances No. 1, 2, and 8 from Group 1, which were clearly assigned to Group 2 already at the first LDA steps. This proves the validity of the mathematical model developed based on physicochemical properties for classifying substances from the class of opioid analgesics into groups.

Thus, in order to determine whether new compounds not included in the research list belong to one of the three groups, it would be necessary to:

It may happen that, when situated relatively far from all the centroids, the determined point will not belong to any of the three groups.

To directly calculate the classification index in order to assign the evaluated compounds to one of the three groups, the CF1-CF3 classification functions can be used to determine which group is most likely to be assigned to a xenobiotic. There are as many classification functions as there are groups:

СF1 = 0.276×MW + 22.406×LogP – 23.622×LogD7.4 + 0.906×Polar – 0.160×TPSA – 66.585

СF2 = 0.137×MW + 20.566×LogP – 21.906× LogD7.4 + 2.657×Polar – 0.105×TPSA – 88.554

СF3 = 0.034×MW + 11.248×LogP – 15.200×LogD7.4 + 3.963×Polar + 0.013×TPSA – 97.835

For each new or known compound, all three functions are calculated, and the compound is assigned to the group for which the classification function is most important.

When new CNS-active substances and xenobiotics are being identified based on the proposed physicochemical criteria, they can be classified into one of three groups: substances with a hazard risk, those with a high hazard risk, and those with an extremely high hazard risk.

The operability of the created mathematical model was assessed using the classification of substances in the training sample. The obtained results of the classification of substances, which differ from the initial classification into groups in the training sample, are shown in Table 10.

According to the results of the conducted classification, a potentially greater hazard was predicted for 9 out of 53 model substances than initially expected. The effectiveness of the proposed model for assessing the potential health hazard of CNS-active substances was demonstrated by the example of propoxyphene. This drug was initially classified as Group 1; however, further analysis clearly allocated it to Group 2, i.e., substances with a high hazard risk. Indeed, since 2009, this drug has been designated as hazardous in a number of countries and discontinued due to the high risk of fatal overdoses and cardiac arrhythmias.2 Table 10 also shows that certain drugs and substances detected in illicit trafficking may carry certain risks when used uncontrolled, which is consistent with the literature data on their toxic effects [22].

Diphenoxylate and buprenorphine, classified as Groups 2 and 3, respectively, can be considered as pharmacological exceptions. Diphenoxylate is used for gastrointestinal disorders, with its action being directed at opioid receptors located in the intestine. Therefore, the excessively high lipophilicity of diphenoxylate prevents its action on the CNS [23]. Buprenorphine, which belongs to the family of highly active Bentley compounds, is a powerful drug that could be just as dangerous as its related etorphine. Nevertheless, pharmacologically, it exhibits the properties of an opioid receptor agonist, posing no significant threat to human life. The results obtained emphasize that the proposed model based on the physicochemical properties of substances cannot take into account the peculiarities of their effect on CNS targets, such as receptors and enzymes. Hence, false positive and false negative results are possible in the resulting solutions. At the same time, the correct classification of most model compounds into groups according to their hazard potential shows the determining effect of the physicochemical properties of narcotic analgesics on their toxic action.

Using the example of representatives of one class of CNS-active substances, we confirmed the relationship between the potential health hazard of narcotic analgesics and their physicochemical properties, which determine the redistribution of such substances from the central bloodstream into the CNS tissues. The molecular weight, lipophilicity, ability to acid-base dissociation, polarity and polarizability of CNS-active substances, as well as their ability to bind to target receptors, equally determine the specific features of their toxic effects. Therefore, these physicochemical properties can serve as criteria for determining the health hazard potential of CNS-active xenobiotics.

The method of linear discriminant analysis with step-by-step inclusion of variables from a training sample of representatives of one class of CNS-active substances (opioid analgesics), divided into three groups according to their hazard level, was used to derive a mathematical model for classifying psychoactive substances based on their physicochemical properties. The model included two dF1 and dF2 discriminant functions as linear combinations of the values of physicochemical descriptors, allowing us to satisfactorily divide the considered groups of narcotic analgesics among themselves according to the criterion of potential threat. The CF1-CF3 classification functions were calculated to assign new compounds to a particular group.

The validity of the proposed mathematical model for assessing the potential hazard of CNS-active xenobiotics was confirmed on the example of classifying a training sample of substances. As a result, 44 substances out of 53 were classified correctly. The classification of seven substances was clarified; two substances proved to be pharmacological exceptions.

The calculated classification functions can be easily embedded as formulas in tabular editors or specialized software and databases for automatic classification of new and known compounds. Provided that the quantitative expression of the power and speed of action of opioid analgesics on the body and other parameters of molecules are determined, it becomes possible to conduct a mathematical analysis in order to establish a quantitative relationship between the characteristics of their biological activity and quantitative structure-activity relationship (QSAR). Of significant interest is the applicability of the results obtained in the study to other classes of CNS-active substances.

Undoubtedly, along with perfection of predictive methods, their results require confirmation by conventional toxicological methods. However, the demonstrated possibility of rapid detection of potentially hazardous xenobiotics will be useful for the prevention of mass poisoning when such CNS-active substances enter illicit traffic.

1. Rembovsky VR, Radilov AS, Dulov SA, Nikolaev AI. Assessment of the degree of exposure to xenobiotics based on the spectral characteristics of water-protein complexes of blood plasma. Methodological recommendations 12.11. Moscow: FMBA; 2012 (In Russ.). 2. Gandey A. Physicians Say Good Riddance to Worst Drug in History.2011. https://www.medscape.com/viewarticle/736718?src=mp&spon=25&form=fpf (Available from: 16 Jun 2024).