Using Delaunay triangulation and Voronoi tessellation to predict the toxicities of binary mixtures containing hormetic compound

Concentration addition (CA) was proposed as a reasonable default approach for the ecological risk assessment of chemical mixtures. However, CA cannot predict the toxicity of mixture at some effect zones if not all components have definite effective concentrations at the given effect, such as some compounds induce hormesis. In this paper, we developed a new method for the toxicity prediction of various types of binary mixtures, an interpolation method based on the Delaunay triangulation (DT) and Voronoi tessellation (VT) as well as the training set of direct equipartition ray design (EquRay) mixtures, simply IDVequ. At first, the EquRay was employed to design the basic concentration compositions of five binary mixture rays. The toxic effects of single components and mixture rays at different times and various concentrations were determined by the time-dependent microplate toxicity analysis. Secondly, the concentration-toxicity data of the pure components and various mixture rays were acted as a training set. The DT triangles and VT polygons were constructed by various vertices of concentrations in the training set. The toxicities of unknown mixtures were predicted by the linear interpolation and natural neighbor interpolation of vertices. The IDVequ successfully predicted the toxicities of various types of binary mixtures.

Delaunay triangulation (DT) is a technique for creating a mesh of contiguous, nonoverlapping triangles from a dataset of points 1 . The DT is the geometric dual of the voronoi tessellation (VT) 2 . A Voronoi diagram splits up a plane based on a set of original points. Each polygon, contains an original point and all areas that are closer to that point than any others 3 . VT and DT have been applied in many areas of mathematics and the natural sciences [4][5][6] . For example, people used linear interpolation (based on DT) and natural neighbor interpolation (based on VT) to predict the climate data or geographic data 7,8 . Natural neighbor interpolation and linear interpolation are the methods of spatial interpolation. Since data acquired in practical production are always limited, spatial interpolation is an effective way to remedy data. Mixture toxicity researching are also faced with the same problem, it is impossible to test the toxicity of all mixtures. This is because multicomponent mixture is a very complex system, the concentration and proportion of components are infinite 9 .
People often use concentration addition (CA) and independent action (IA) to predict the toxicity of the mixture [10][11][12] . In particular, CA was proposed as a reasonable default approach for ecological risk assessment of chemical mixtures 13 . CA can be expressed as shown in equation (1): where n refers to the number of components of mixture, c i denotes the concentration of the ith component in the mixture that elicits x% effect, and EC x,i is the effective concentration of the component i that provokes the same effect (x%) when applied alone. The fraction c i /EC x,i is often termed a 'toxic unit (TU i )' , and CA is hence also known as 'toxic unit summation' 14 .
Scientific RepoRts | 7:43473 | DOI: 10.1038/srep43473 CA is based on the premise that all components have similar mechanisms of action (MOA) 15 . IA assumes that components of mixture act dissimilarly 16 . However, toxicologists often know little about the MOAs of substances 17 . There is an obvious defect of CA, the only condition for CA predicting the effective concentrations (EC x ) is that all components have definite effective concentrations at the given effect (x%). Therefore, CA cannot predict the toxicity of mixture at some effect zones (predictive blind zones) when the mixture and the components have different effects 18 , such as some compounds induce hormesis 19 . Hormesis is a concentration-response phenomenon that is characterized by low-dose stimulation and high-dose inhibition 20 . IA can predict mixture effect when all or some of the chemicals induce hormesis in theory, but in practice, some people think IA will lose its conceptual framework if negative values for single chemical effect were involved 21 . The calculation for the IA model can be performed according to equation (2): where E(c mix ) is the predictive effect of a mixture with a total concentration of c mix , c i is the individual concentration of ith compound in the mixture and E(c i ) is the effect of this concentration if the compound is applied singly. In fact, it isn't uncommon for the toxicity of a mixture deviates from such predictions when antagonism or synergism occurs. If the toxicological interaction, antagonism or synergism, is identified by CA at a certain mixture ratio or concentration level, CA can not certainly predict the toxicity accurately at other mixture ratio or concentration level. Because the interaction is mixture ratio-dependent and concentration level-dependent 22,23 . For example, we have determined the mixture toxicities of 1-ethyl-3-methylimidazolium chloride ([emim]Cl) and metalaxyl (MET) on Vibrio qinghaiensis sp.-Q67 (V. qinghaiensis), and found that the mixture exhibited synergism at one mixture ratio, but showed antagonism at another ratio 24 . It means that toxicological interaction of chemicals does occur neither uniformly nor predictably. Besides CA and IA, some researchers tried to predict the toxicity of mixture through using the QSAR-based approaches to predict the toxicities of single compound and then compute the mixture descriptors [25][26][27] . However, QSAR is very successful in dealing with individual compounds, but how to calculate rationally the mixture descriptors is still in dispute.
In this paper, we developed a new method for the toxicity prediction of various types of binary mixtures ( Fig. 1), an interpolation method based on the DT and VT as well as the training set of direct equipartition ray design (EquRay) mixtures, simply IDV equ . In order to obtain the required data, time-dependent toxicity of two ILs, three pesticides and binary pesticide-IL mixtures on V. qinghaiensis were determined. Because we aim to test the predictive capability of the methods when toxicological interaction occurs, some results of our previous researches were used. They are the mixture toxicities of [emim]Cl and MET, and mixture toxicities of 1-ethyl-3-methylimidazolium bromine ([emim]Br) and MET 24 .

Results
Concentration-response curves (CRCs) of five single substances. The concentration-response curves (CRCs) of three pesticides and two ILs at seven exposure times were displayed in Fig. 2 70 and characteristic parameters (zero effect point (ZEP), minimum effective concentration (EC min ) and E min for J-shape CRCs) 28,29 of pesticides and ILs at seven exposure times were reported in the Table S1. The experimental concentration-inhibition points, fitted CRCs of five compounds and 95% confidence intervals (CI) at seven exposure times were shown in Figure S1.  Table 1. For example, the R3 ray in [epy]Br-DOD systems, has no interaction at the exposure times of shorter than 4 h and the effect levels of larger than 60%, the corresponding ERR values basically close to zero and exihibit antagonism at the times of longer than 4 h and the effects of less than 60%, the corresponding ERRs significantly less than zero.
Performance of IDV equ . The accuracy rate (AR), R 2 and RMSE of IDV equ predicting the toxicities of various mixtures in eight binary mixture systems by the leave-many-out cross-validation (LMOCV) and leave-one-out cross-validation (LOOCV) analysis were displayed in Table 2 In Table 2  plots by LOOCV were showed in Figure S6a and S6b. The predictive ability of LinIP and NeiIP is almost the same. However, the results from LMOCV and LOOCV are different. The ARs from LOOCV analysis are, as a whole, lower than those from LMOCV.

Discussion
IDV equ is independent of the chemical MOA. When CA or IA is used to predict the toxicity of the mixture, people need to know the MOA of various mixture components. However, the MOA of most chemicals still remain unknown, it will result in difficulties to choose CA or IA as the reference models 30 . Furthermore, realworld mixtures are made up of chemicals with both similar and dissimilar MOA 31 . Although the IDV equ requires some data of mixture toxicity, it doesn't require the MOA of each chemical. IDV equ is based on the experiment datas, it is more reliable because no model can guarantee the correctness of the prediction.    32,33 and contaminants always co-occour in ecosystems 16 . The results indicated that IDV equ can predict the toxicity of binary mixtures which contain hormetic chemical well. The results of this paper also show that, even if interaction occurs, IDV equ can also well predict the toxicity at other ratios or concentrations. But if interaction is identified by CA at a certain ratio or concentration, whether CA can predict the toxicity accurately at other ratios or concentrations is unknown. Our previous research 24 and the results of this paper showed that toxicological interaction in the pesticide-IL mixture system is ratio-dependent and time-dependent. An empirical model was formulated by Jonker et al. 34 to determine the magnitude of interaction, it requires experimental designing which covers all ratios and concentration for optimal performance. Unforunately, the combination ratio and concentration of mixtures is infinite 16,35,36 . In this paper, IDV equ predict the interaction just need the toxicity of five mixture rays and two single compounds. It will minimising the experimental investigations and the consequent consumption of time and resources.
Largest interactions does not always occur at equitoxic ratios. Some people think the largest interactions most often occur at equitoxic ratios. But the results of this paper are not consistent with this hypothesis. In [epy]Br-DOD system, the maximum ERR occur at R3 (0.25 h), R2 (2 h), R3 (4 h), R3 (6 h), R4 (8 h), R4 (10 h) and R4 (12 h). In [epy]Cl-DOD system, the maximum ERR occur at R3 (0.25 h), R3 (2 h), R4 (4 h), R4 (6 h), R3 (8 h), R4 (10 h) and R4 (12 h). According to EquRay, the mixture ratio of R3 is 1:1 base on EC 50 (equitoxic ratio), but the maximum ERR don't always occur at R3. There are a lot of researches showed that the largest interactions don't often occur at equitoxic ratios. For example, people tested the binary mixture of zinc and copper on Tympanotonus fuscatus and found that the mixture (1:4 mixture ratio base EC 50 ) exhibited strong antagonism, but the mixture showed additive at equitoxic ratio 37 . Cedergreen 38 invertigated the binary mixture of mecoprop and terbuthylazine on the floating plant Lemna minor and found that antagonism was larger for mixtures with higher proportions of mecoprop.
There are no predictive blind zones for IDV equ . CA cannot predict the toxicity of mixture at some effect zones (predictive blind zones) when the compounds share a different range of effects. The predictive blind zones contains a variety of types: partial mixture components can cause hormesis; all mixture components can horemsis but their E min is not equal; the maximum effect of a mixture is different from that of a single component. In these cases, IDV equ is able to predict the toxicity of the mixture at all effect zones. Although a generalized concentration addition (GCA) model can predict the toxicity of mixtures when some chemicals have a smaller maximum effect level than others 39 , it need assume that all CRCs were fitted by Hill functions and slope parameter was one. This limits the use of the GCA model.
People always hope that a model can be applied to all situations. For example, some people proposed a modified model and automated fitting procedure to describe CRCs with multiphasic features 40 . The application of IDV equ is not restricted by the CRC of single compounds and mixtures. There are no predictive blind zones for IDV equ , it should be a universally applicable method for predicting the mixture toxicity. Theoretically, LinIP and NeiIP method can predict the toxicity of multi-component mixtures. In fact, we have studied the toxicity of multi-component mixtures in our previous works. For example, we used UD-Ray method to select the representative mixtures from a lot of mixtures rays in the multi-component mixture system 41 . But the interpolation methods to predict the toxicity of multi-component mixtures need to be further validated. We expect IDV equ to be a tool useful for experimentalists and analysts interested in the study of mixture toxicity, because people are rarely exposed to a single hazardous chemical 42 . Cl-MET and [epy]Cl-SIM. EquRay was employed to design the basic concentration compositions of five binary mixture rays 43 . The mixture ratios (p i,j ) 44 , the ratio of the concentration of the jth component in the ith ray to the total concentration of the ray, of various components in 30 mixture rays and the concentrations of stocks were listed in Table 3.

Materials and Methods
Determine the toxicities of single components and mixture rays. A time-dependent microplate toxicity analysis was given in our previous works 45 . The toxic effects of single components and mixture rays at different times and various concentrations were determined by it. The setup of controls and treat-group in the 96-well microplate is designed according to Fig. 1 in the literature 46 . The 12 concentration gradients were calculated according to equation 3 in in our previous work 47 . The freeze-dried V. qinghaiensis was purchased from Beijing Hamamatsu Corp., Ltd. (Beijing, China). The relative light unit (RLU) of every well was determined on the Power-Ware microplate spectrophotometer (American BIO-TEK Company) at 22 ± 1 °C. During exposure, readings were taken at 0.25, 2, 4, 6, 8, 10 and 12 h. Inhibition ratio of bioluminescence was used to characterize the toxicity, noted as E: where I 0 is the average RLU of V. qinghaiensis exposed to the experimental group and I indicates the average RLU of V. qinghaiensis exposed to controls. 1, 2, …, 30; j = 1, 2, 3, 4, 5)    Using LinIP and NeiIP to predict the toxicity of unknown mixtures. The toxicities of unknown mixtures were predicted by the LinIP and NeiIP of vertices. The sketch-map of LinIP and natural NeiIP for predicting the toxicity of binary mixture is displayed in Fig. 4. Firstly of all, define the point P like this: (X i ,Y i ,Z i ), X i means concentration of one substance, Y i means concentration of another chemical, Z i means the toxicity of the binary mixture. LinIP is used to interpolate the toxicity of query point. The operation steps of LinIP are as follows 7 : each triangle covers an area and the toxicity of query points in this area is predicted based on the triangle that covers it. Let P 1 , P 2 , and P 3 be the three vertices of the triangle, located at P 1 = (X 1 ,Y 1 ,Z 1 ), P 2 = (X 2 ,Y 2 ,Z 2 ) and P 3 = (X 3 ,Y 3 ,Z 3 ). Suppose that the plane which is defined by the three points is given by

Mixture ray
We will obtain a linear system of equations by inserting the X, Y and Z values of each of the three known points into this equation 50  where the coefficients a, b and c of the equation (4) can be found when we solve this system of equations. The value of Z for any query point within this triangle will be available after the above steps. Although voronoi tesselation of given scattered point set has been constructed first, when a query point is inserted a new voronoi polygon is created which overlaps the original voronoi 48 . The intersection of the new Voronoi polygon with the original Voronoi diagram is defined 51 as follows: i i  where V i is the Voronoi polygon of the ith natural neighbor x i . It is said that the new Voronoi polygon V(q) overlaps some area from the neighboring Voronoi tiles V i of the query point. Then, the weight of ith natural neighbor, W i , is defined by is the area of each intersection Vi(q). The Area [V(q)] is the total area of the polygon V(q). Then, the NeiIP of the query point q and is defined by i i i where f(x i ) is the function value in the ith natural neighbor point. f(q) is the function value of the query point.
Validation of the predictive capabilities of IDV equ . LMOCV and LOOCV were selected to determine the performance of IDV equ . LMOCV randomly splits the dataset into training and validation data, approximately 2/3 of the data was allocated to the training set while the remaining 1/3 was allocated to the test set 52 . In this paper, mixture toxicity (n = 60) were randomly divided into two disjoint parts; 1/3 of the combined toxicity (n = 20) was allocated to the test set, 2/3 of the combined toxicity (n = 40) and single toxicity (n = 24) were combined resulting in a training set. The holdout method was repeated 20 times with randomly selected training and holdout sets 53 . LOOCV is a particular case of LMOCV where M = O (one). In this paper, the toxicity of a mixture ray from the mixture system is selected as a validation set. The toxicity of remaining four mixture rays and two individual substances are used as the training set. The purpose of this is to test whether the two methods can accurately predict the complete CRCs.
The predicted values were compared with 95% CI of observed value. If CI overlaps predictions, we believe that the predicted value is correct.The root mean square error (RMSE), coefficient of determination (R 2 ) and accuracy rate (AR) were used to assess the performance. Accuracy rate (AR) is the number of query point be predicted correctly divided by the number of total query point.

CRCs fitting and toxicity interaction characterization. Monotonic S-shaped CRCs were fitted by Logit
or Weibull function 54 . A five parameters logistic equation (equation 8) was chosen to describe the non-monotonic J-shaped CRCs 55 . CA was used to predicted the toxicity of mixtures rays. where E min refers to the minimum effect or maximum stimulatory effect, ε dn refers to the concentration at the effect of E min /2 in the falling section (negative slope), β dn refers to the slope at the point (ε dn , E min /2), ε up refers to the median effective concentration, β up to the slope at the point (ε up , 50), EC x refers to concentration. ERR x (equation 9) was used to quantify the toxicological interaction (synergism or antagonism) in mixtures 24 . Considering the CI, the value of ERR x at a specific effect (x) can be computed as follows: x CI prd prd where E CI is the effect (toxic response) corresponding to the upper limit (when antagonism occurs) or lower limit (when synergism occurs) of CI, E prd is the effect value predicted by an additive reference model such as concentration addition (CA) at the same concentration (EC x ). When ERR x > , = , and < 0, say the mixture at the effect of x being synerigism, additive action, and antagonism.