Mitogen-activated protein kinase (MAPK) modulates extracellular signals to control several cell functions (from proliferation and survival to differentiation and senescence)^[1,2]. One of the most studied MAPK pathways is the extracellular signal-regulated kinase (ERK) pathway. ERK is a subgroup of MAPKs that is activated by external factors such as growth factors and mitogens^[3]. The cascade of ERK signalling is triggered by the activation of receptor tyrosine kinases and flows through RAS GTPase^[4]. RAF kinases are key components of this pathway: once RAS is turned on, it recruits and activates proteins necessary for the propagation of growth factor and other receptor signals, such as RAF. There are three RAF proteins in mammals, ARAF, BRAF, and CRAF, and they can all activate MAP kinase (MEK) just upstream of ERK^[5]. The over-expression or mutations of the components of the RAS–RAF–MEK–ERK–MAP kinase pathway has been found to occur in up to 30 % of human cancers^[6]. In different cancer types, the mutation of RAS or RAF causes the hyper activation of ERK followed by unrepresents a potential approach for the treatment of cancer^[8]. In particular, BRAF somatic missense mutations have been found to occur in 66% of malignant melanomas and at lower frequency in a wide range of human cancers^[5,9].

  Among the BRAF mutations, a point mutation resulting in substitution of glutamic acid for valine, also known as BRAFV600E, is the most common change^[6]. In the recent past, a number of small molecules inhibiting mutant RAF kinases have been discovered, and some of these have now advanced into clinical trials. For example, Sorafenib^[10] (Figure 1) is the most studied inhibitor for this class of kinases. It was approved by FDA in 2005 for the treatment of renal cell carcinoma and in 2007 for the treatment of hepatocellular carcinoma, and is still undergoing multiple clinical trials in other types of cancer^[11,12].

Figure 1: Structures of Sorafenib (Bay 43-9006; Nexavar).

  The experimental measurement of the inhibition activity of chemicals is difficult, expensive and time-consuming. QSAR based analysis can be used as a tool to screen or filter anti-cancer drug candidates, before they are subjected to more intensive calculations, such as docking, or to experimental in vitro measurement of activity and finally under in vivo conditions^[13]. In this study, a model was built to predict the half maximal inhibitory concentration (IC₅₀) of ^V600EBRAF-dependent ERK phosphorylation in WM266.4 cells. This model made use of descriptors based on simplified molecular input-line entry system (SMILES)^[14-17], calculated with the CORAL software. Then, chemical information obtained from the calculation of these descriptors was used for drawing up hypothesis of molecular design of Sorafenib derivatives.

  The aim was to obtain molecules able to inhibit 50% of ERK phosphorylation in WM266.4 melanoma cells at lower concentrations compared to the template molecule.

Methods

Data
  The dataset of 142 chemical structures used in this study was taken from the literature^[18,19]. All these compounds shared a common scaffold (Figure 2). The dataset contained chemical structures represented as SMILES strings, and experimental values relative to the inhibition of the phosphorylation of extracellular signal-regulated kinases (p-ERK).

Figure 2: General structure of target series.

  The selected endpoint was IC50 of V600EBRAF-dependent^[20] ERK phosphorylation in WM266.4 melanoma cells (IC₅₀, pERK)^[1,14]. IC₅₀ values were expressed as negative decimal logarithm (pIC₅₀). In this study ACD ⁄ ChemSketch^[21], was used to build up SMILES of each structure.

CORAL and Optimal descriptors
  CORAL is freely available standalone application software for building up regression or classification QSAR models based on the Monte Carlo technique^[22]. The calculations are repeated several times for various splits into training, calibration and external validation sets. Optimal descriptors for constructing QSAR models are based on graph or SMILES. The molecular graph includes hydrogen suppressed graph (HSG), hydrogen filled graph (HFG), and graph of atomic orbital’s (GAO).

  The SMILES-based optimal descriptors include a group of local and one of global attributes. Sk, SSk, and SSSk are local descriptors represented by a sequence of atoms and bonds present in the SMILES string^[23].

The SMILES “ClCCCBr”, Sk, SSk, and SSSk can be represented as follows:
ClCCCBr → Cl + C + C + C + Br ( Sk )
ClCCCBr → ClC + CC + CC + CBr( SSk )
ClCCCBr → ClCC +CCC+ CCBr( SSSk )

BOND, PAIR, NOSP and HALO are global SMILES attributes, involving the presence of the following specific items in the target SMILES:
BOND refers to the presence/absence of double (=), triple (#), and stereo chemical (@) bonds;
PAIR indicates the co-occurrence of two elements among the followings: F, Cl, Br, I, N, O, S, P, #, =, @;
NOSP indicates the presence/absence of nitrogen, oxygen, sulfur, and phosphorus;
HALO refers to the presence/absence of halogens.

  A more detailed explanation of CORAL descriptors has been provided by^[24]. All the SMILES-based local descriptors and BOND, NOSP and HALO global attributes were selected to build the regression models.

Table 1: The statistical characteristics of QSAR model of pIC₅₀ p-ERK inhibitors.

  Data were randomly split three times into training (60%), calibration (20%) and validation sets (21%). Table 1 The training and calibration sets were used for building the QSAR model and the validation set was used as external set for estimating the predictive potential of the developed model^[25]. In order to develop a model with a good predictive potential, the preferable parameters of the Monte Carlo optimization, the threshold (T*) and the number of epochs (N*) that give the best statistics for the calibration set should be defined.

  The threshold is a tool for classifying codes as both rare (and thus likely less reliable features, probably introducing noise into the model) and common features, which are used by the model and labelled as active. The optimal descriptors are calculated with the correlation weights (CW) only of active features, excluding those related to rare ones.

  The Nepoch is the number of cycles (sequence of modifications of correlation weights for all codes involved in model development) for the optimization^[26]. For modelling activity of potential p-ERK inhibitors the descriptor of correlation weights (DCW) was calculated as follows:
DCW (Threshold, Nepoch) = Σ CW (Sk) + Σ CW (SSk) + Σ CW (SSSk) + Σ CW (BOND) + Σ CW (HALO) + Σ CW (NOSP)    (1)
where CW(SAk) is correlation weight for a molecular or structural feature (SAk) extracted from k-th SMILES;

  The endpoint pIC₅₀ is a function of DCW according to the following equation:
Endpoint = C0 + C1 x DCW (T, Nepoch)    (2)
Where C0 and C1 are the intercept and the slope for the training and calibration set, respectively.

  The role of molecular features or structural attributes (SAk) extracted by SMILES can be defined by the value of CWs: SAk with positive CW are promoters of pIC₅₀; SAk with negative CWs are responsible for the endpoint decrease; if there are both positive and negative values of CW (Sk), then that SAk has an undefined role. The applicability domain of the models can be defined according to the distribution of the molecular attributes extracted from the k-th SMILES.

  A SMILES falls into the applicability domain if the following condition occurs:
Defect SMILES < 2* (Defect SMILES)̅    (3)
Where Defect SMILES is the measure of the statistical (probabilistic) quality of molecular features extracted from k-th SMILES, and (Defect SMILES)̅ is the average of these values for the training and calibration sets. Defect SMILES is calculated as follows:
Defect SMILES = Σ≣ {Defect({SA}_k)}    (4)
Defect ({SA}_k) = |P_(TRN(SAk)-) P_TST(SAk) |/(N_(TRN(- SAk)+) N_TST(SAk) )     (5)

  Where PTRN (SAk) is the probability of the presence of SAk in the SMILES of the training set and PTST (SAk) is the probability of SAk in the SMILES of the calibration set. NTRN (SAk) is the number (frequency) of SMILES containing SAk in the training set and NTST (SAk) is the number of SMILES containing SAk in the calibration set. The ideal situation occurs when the probabilities of SAk are the same in the training and in the calibration sets (Defect (SAk) is equal to zero). Conversely, if SAk is absent in the calibration set, the Defect (SAk) is maximal. Thus, the results from equation (4) can be used to classify the active attributes, and the Defect SMILES defines the domain of applicability for SMILES^[26].

Results and Discussion

The models for three random splits are the following:
Split 01 pIC₅₀ = 1.4942 (± 0.0237) + 0.0379 (± 0.0001) * DCW (3, 3) (07)
Split 02 pIC₅₀ = -1.0822 (± 0.0222) + 0.0342 (± 0.0001) * DCW (5, 4) (08)
Split 03 pIC₅₀ = -1.8219 (± 0.0236) + 0.05103 (± 0.0002) * DCW (2, 6) (09)

Figure 3(a), Figure 3(b), Figure 3(c) shown a graphs plotting of the calculated versus experimental pIC₅₀ values.

Figure 3a: Graphical representation of built QSAR model for split 01.

Figure 3b: Graphical representation of built QSAR model for split 02.

Figure 3c: Graphical representation of built QSAR model for split 03.

  The values of the parameters of statistical analysis of p-ERK inhibition models built on the three random distributions of data into training, calibration, and validation sets are collected in Table 2. n is the number of compounds in each set; r² is the coefficient of determination; q² is cross validated r²; MAE is mean absolute error; s is the root-mean-square error; F is the Fischer F ratio; T and N are the preferable values for Threshold and Number of epochs, respectively. CR_p^2 indicates if the models developed are obtained by chance or not, based on the Y-randomization test. For an acceptable model, CR_p^2 should be greater than 0.5^[27,28]. r² range from 0.78 to 0.92, and q² range from 0.75 to 0.92 (for training, calibration and test sets). According to the criteria defined by Tropsha^[26-29], a model has high predictive power if the following conditions are fulfilled:
q² > 0.5
r² > 0.6
(r² - r²0)/r² < 0.1 or (r² - r′0²)/r² < 0.1
0.85 ≤ k ≤ 1.15 or 0.85 ≤ k′ ≤ 1.15
|r0²−r′0²| < 0.3

  Where r0² and r′0² are squared correlation coefficients for regression through the origin, calculated between predicted versus experimental values and between experimental versus predicted values, k and k′ are the slopes in the former and later cases respectively; q² is calculated for the training sets, while all other criteria are calculated for the validation sets^[26]. Table 3 contains the values for the criteria iii, iv, and v. Y-randomization test for all models showed that these are not chance correlations, since the CR_p^2 is larger than 0.5 (Table 2). All the models successfully fulfill the criteria proposed by Tropsha for predictive ability. 86 out of 142 dataset’s compounds assigned to the training set are used to develop the QSAR models, 28 compounds are included in the calibration, and 28 compounds, not involved in the model development, are kept out for the external validation. Based on r² and q² values, the three models of equations (07), (08) and (09) have similar statistical behaviours. Indicating that the random selection does not affect the results. The predictability of the models was also checked calculating the statistical parameters r_m², average r_m², and Δr_m² on the calibration sets (Table 3), as proposed in the literature^[29,30]. According to the literature, a model has predictive potential if r_m² > 0.5, average r_m² > 0.5, and Δr_m² < 0.2^[31]. The presented results reveal that the predict ability of all models is very good, especially for the model of equation (09), with q² and r_m² values equal to e 0.92 and 0.88, respectively (Table 3).

Table 2: The numerical values of (r²-r′0²)/r²; (r²- r0²)/r²; k; k′ and |r0² – r′0²|.

Table 3: The numerical values of r_m²; (r_m²)̅; Δr_m².

  The list of all SAk, with the correlation weights for the three splits from the Monte Carlo optimization process of the built QSAR model is given in Table S2, S3 and S4. Based on the best model in validation obtained from split 2, some SMILES based descriptors have positive influence on pIC₅₀ and therefore cause an increase of pIC₅₀ value^[32-35]. : =........... (Molecule containing double bond), N........... (Molecule contains nitrogen atom), O........... (Molecule containing oxygen atom), c...O....... (Molecule contains oxygen atom bonded to aromatic carbon atom), c...N....... (Molecule containing nitrogen atom bonded to aromatic carbon atom), c...S....... (Molecule containing sulfur atom bonded to aromatic carbon atom), Cl.......... (Molecule containing chlorine atom), BOND10000000: (molecule contains fluorine atom) etc.

  The mechanistic interpretation obtained from the analysis of SAk can be used in the search and computer aided design of novel p-ERK inhibitors candidates with desired pIC₅₀ values. Figure 4 shows the structures of proposed kinase inhibitors derivatives, the structure Sorafenib was used as a template for molecular design by added the SMILES notation c...S......., BOND10000000 and c...O....... defined as promoters of increase. The molecule E has a higher value of pIC₅₀ in comparison to other molecules.

Figure 4: The molecular design of possible kinase inhibitors using the QSAR model of split 2 and SAks calculated using the Monte Carlo method

Table 4: SMILES notation and pIC₅₀ values calculated by using model of split 2 for kinase inhibitors derivatives designed by using the results of QSAR model obtained in this study.

  Supplementary materials section contains technical details of the models for splits examined in this work.

Conclusion

  In this work three QSAR models were developed to predict ^V600EBRAF-dependent ERK phosphorylation of 142 molecules as kinase inhibitors, using SMILES based optimal descriptors. The models showed acceptable predictive capability on three random split of data. The analysis of the structural features (obtained from SMILES) with their positive and negative correlation weights allowed designing possible pERK inhibitors with an increased activity compared to Surafenib.

  The activity modulating effect of the structural features may serve to draw up preliminary hypothesis of novel anticancer drug candidates, and also to provide a better understanding the drugs mechanisms of action.

Acknowledge:
  This research has been partially done during a research visit of the Corresponding author to Politecnico di Milano and Mario Negri Institute in Italy, in the PROSIL Life + research project of the European Union (LIFE12 ENV/IT/000154).

Conflict of interest:
  The authors have no conflict of interest to declare.

Supplementary material

Table S1. Experimental and calculated of pIC₅₀ p-ERK inhibitors for three splits; training (Tr); calibration (C); and validation (V) sets.

Table S2. Correlation weights for DCW calculation for split 1; obtained in three probes of the Monte Carlo optimization.

Table S3. Correlation weights for DCW calculation for split 2; obtained in three probes of the Monte Carlo optimization.

Table S4. Correlation weights for DCW calculation for split 3; obtained in three probes of the Monte Carlo optimization.

References