Head to head evaluation of second generation ALK inhibitors brigatinib and alectinib as first-line treatment for ALK+ NSCLC using an in silico systems biology-based approach

Around 3–7% of patients with non-small cell lung cancer (NSCLC), which represent 85% of diagnosed lung cancers, have a rearrangement in the ALK gene that produces an abnormal activity of the ALK protein cell signaling pathway. The developed ALK tyrosine kinase inhibitors (TKIs), such as crizotinib, ceritinib, alectinib, brigatinib and lorlatinb present good performance treating ALK+ NSCLC, although all patients invariably develop resistance due to ALK secondary mutations or bypass mechanisms. In the present study, we compare the potential differences between brigatinib and alectinib’s mechanisms of action as first-line treatment for ALK+ NSCLC in a systems biology-based in silico setting. Therapeutic performance mapping system (TPMS) technology was used to characterize the mechanisms of action of brigatinib and alectinib and the impact of potential resistances and drug interferences with concomitant treatments. The analyses indicate that brigatinib and alectinib affect cell growth, apoptosis and immune evasion through ALK inhibition. However, brigatinib seems to achieve a more diverse downstream effect due to a broader cancer-related kinase target spectrum. Brigatinib also shows a robust effect over invasiveness and central nervous system metastasis-related mechanisms, whereas alectinib seems to have a greater impact on the immune evasion mechanism. Based on this in silico head to head study, we conclude that brigatinib shows a predicted efficacy similar to alectinib and could be a good candidate in a first-line setting against ALK+ NSCLC. Future investigation involving clinical studies will be needed to confirm these findings. These in silico systems biology-based models could be applied for exploring other unanswered questions.

Systems biology-based models were created using the Therapeutic Performance Mapping System (TPMS) to investigate the molecular Mechanisms of Action (MoA) of brigatinib and alectinib towards the modulation of ALK+ NSCLC. TPMS is a validated top-down systems biology approach that integrates all available biological, pharmacological and medical knowledge by means of pattern recognition models and artificial intelligence to create mathematical models that simulate in silico the behavior of human physiology. The methodology employed and detailed herein has been previously described [60] and applied elsewhere [59,115,116,121]. Biological maps were transformed into a mathematical model capable of both reproducing existing knowledge and predicting new data. TPMS technology uses a set of artificial intelligence algorithms to generate the human physiology over the human biological network [121][122][123][124].

Human protein network (HPN) and truth table construction
A human protein network (HPN) was created in order to obtain the MoAs by including information from many public and private databases (KEGG [125,126], REACTOME, [127] INTACT, [128] BIOGRID, [129] HPRD, [130] and TRRUST [131] and information extracted from scientific literature. For the construction of the truth table, a selected collection of known input-output physiological signals considered the "truths" were collated into a table  (Supplementary Table 7) and was used for training the models [132]. The truth table was based on a compendium of different databases that contain biological and clinical data [133,134] and provides biological and pharmacological input-output relationships (such as drugindication pairs). Information relating biological processes (adverse drug reactions, indications, diseases and molecular pathways) to their molecular effectors, i.e., each one of the proteins involved in the physiological process, was extracted from the biological effectors database (BED) (Anaxomics Biotech SL, [60]). The biological or pathological conditions under study were also included in the truth table and molecularly characterized through specific scientific literature search and hand-curated assignment of proteins to the conditions (Supplementary Table 2). The obtained final models had to be able to reproduce every rule contained in the truth table, and we defined the error of a model as the percentage of all the rules with which the model does not comply, while the accuracy was defined as the percentage of all the rules complied with.

Modelling strategies
Two complementary modelling strategies were used, (a) TPMS Artificial Neural Networks (ANNs) [59] and (b) TPMS Sampling-based Methods [60], to compare the efficacy of the drugs (defined as their targets, see Supplementary Table 3) and to compute the MoA models.
(a) ANNs are supervised algorithms that identify relations between proteins (e.g., drug targets) and clinical elements of a protein network [59,118,120,135,136] by inferring the probability of the existence of a specific relationship between two or more protein sets, based on the validation of the predictive capacity of the model towards the truth table. The learning methodology used consisted in an architecture of stratified ensembles of neural networks as a model, trained with a gradient descent algorithm to approximate the values of the given truth table. The neural network model used consisted in a Multilayer Perceptron (MLP) neural network classifier. MLP gradient descent training depends on randomization initialization and to avoid random errors 1000 MLPs are trained with the training subset and the best 100 MLPs are used. In order to correctly predict the effect of a drug independently of the number of targets, a different ensemble of neural networks are trained for a different subset of drugs according to their number of targets (drugs with 1 target, 2 targets, 3 targets). Then, the predictions for a query drug are calculated by all the ensembles, and pondered according to the number of targets of the query drug (the difference between the number of targets of the query and the number of targets of the drugs used to calculate each ensemble is used to ponder the result of each ensemble). A cross-validation with the truth table information showed that the accuracy of the described ANNs to reproduce the indications compiled in DrugBank [112,133] is 81.7% for those drugs with all targets in the human biological network.
(b) Sampling-based methods generate models similar to a MLP over the previously constructed HPN, where neurons are the proteins and the edges of the network are used to transfer the information (Supplementary Figure 1) [60]. This methodology was used for describing with high capability all plausible relationships between an input (or stimulus) and an output (or response). Sampling-based methods use optimization algorithms [123] to solve each parameter of the equation, i.e. the weights associated to the links between the nodes in the human protein network. In this approach, the network is limited by considering only interactions that connect drug targets with protein effectors in a maximum of three steps. The values of activation (+1) and inactivation (-1) of the protein targets of the drugs in the truth table were considered as input signals whereas the output is defined as the values of activation and inactivation of the proteins describing the phenotype (as retrieved from the BED). Each node of the protein network receives as input the output of the connected nodes in the direction flow from targets to effectors, weighted by each link weight (Supplementary Figure 1). The sum of inputs is transformed by a hyperbolic tangent function to generate the score of the node (neuron), which becomes the 'output signal' of the current node towards the nodes. The weight parameters are obtained by Stochastic Optimization Method based on Simulated Annealing [123], which uses probabilistic measures derived from the biological evidence to adjust network interaction types and strengths. Since the number of entries in the truth table is always smaller than the number of parameters (link weights) required by the algorithm, any process modelled by TPMS considers a population of different solutions.

Mechanisms of action elucidation
The MoAs obtained with the TPMS simulates potential interactions between drug targets and protein effectors associated to prototype-ALK+ NSCLC patients. In order to validate this approach, the intensity of the model's response, divided in TSignal and number of protein effectors activated, was used to understand the relationships between all potential mechanisms and compare sets of MoAs from different views (Supplementary Figure 1) [60].

Intensity of the response
We defined the "intensity" of the response as follows: 1) the quantity of protein effectors (#) that reach an expected signal sign; and 2) the strength or amount of the output signal reaching the effectors (i.e., a global measure of the output signal, named TSignal).
Given a protein effector "i", which reaches a signal value y i , and v i being the effector sign according to the BED (active or inactive) and n is the total number of effectors described for a phenotype, it was determined:

Number of effectors achieving the expected sign
Assuming that a drug may be able to activate/ inactivate protein effectors reverting a disease/indication model phenotype. Using Dirac's δ (i.e. δ(0) = 1, and zero otherwise), the equation to calculate number of effectors achieving the expected sign for drug indications was defined as:

TSignal
The average output values of the protein effectors. For each effector, it was counted as positive signal if the sign is correct, and negative otherwise. When a drug affects a disease phenotype, v i and y i have opposite sign and it is necessary to change the sign in the corresponding equation:

Sobol sensitivity analysis
An adapted methodology for ensembles of high dimensional algorithms was applied following the definition of Sobol Sensitive Analysis [117]. According to the Sobol terminology, TPMS models can be redefined as follows: and 3B] Where x i is each of the parameters used in the TPMS models. Then, the variation of TSignal for each x i parameter can be expressed as: Consequently, the variation of the simultaneous parameters x i and x j can be estimated as: Using the previous equation descriptions, we measured the impact of varying random parameters over output TSignal in two different approaches, those being local analysis and global analysis [117].

Local sensitivity analysis
Local sensitivity analysis evaluates changes in the model outputs (TSignal) with respect to variations in a single model parameter. This effect was measured in the TPMS-models for both alectinib and brigatinib MoA models (Supplementary Figure 2).

Global sensitivity analysis
In the global sensitivity analysis, all parameters are varied simultaneously over the entire parameter space to measure the effects of their interactions on the model output. Given the high dimensionality of the TPMS ensemble models, this measure has been estimated by a MonteCarlo experiment to introduce random values (noise) in sets of 1200 candidate parameters. These final TSignal effects were measured by altering combinations of parameters in subsets of 1, 2, 3, 4, 5, 10, 15, 20 and 30 parameters simultaneously, from the candidate parameters list (Supplementary Figure 3).

Sensitivity results
Although TPMS-models have about 5000 parameters, only a small percentage of them showed a real impact on the output, which was less notorious in brigatinib than alectinib ( Supplementary Figures 2 and 3). Nevertheless, the impact of some of the protein parameters are of great importance, meaning that TPMS models had to carefully adjust to all the restrictions defines in the truth table, while completing the drug-pathology model. We can see this as most protein parameters are actually part of the ALK+ NSCLC effectors (like P27361, P28482 and P414921, among others), which will definitely have a huge effect on the final TSignal according to its definition in Equation 2.