Tuning the Biological Activity of PI3Kδ Inhibitor by the Introduction of a Fluorine Atom Using the Computational Workflow

As a member of the class I PI3K family, phosphoinositide 3-kinase δ (PI3Kδ) is an important signaling biomolecule that controls immune cell differentiation, proliferation, migration, and survival. It also represents a potential and promising therapeutic approach for the management of numerous inflammatory and autoimmune diseases. We designed and assessed the biological activity of new fluorinated analogues of CPL302415, taking into account the therapeutic potential of our selective PI3K inhibitor and fluorine introduction as one of the most frequently used modifications of a lead compound to further improve its biological activity. In this paper, we compare and evaluate the accuracy of our previously described and validated in silico workflow with that of the standard (rigid) molecular docking approach. The findings demonstrated that a properly fitted catalytic (binding) pocket for our chemical cores at the induced-fit docking (IFD) and molecular dynamics (MD) stages, along with QM-derived atomic charges, can be used for activity prediction to better distinguish between active and inactive molecules. Moreover, the standard approach seems to be insufficient to score the halogenated derivatives due to the fixed atomic charges, which do not consider the response and indictive effects caused by fluorine. The proposed computational workflow provides a computational tool for the rational design of novel halogenated drugs.


Introduction
The inhibition of phosphoinositide 3-kinase (PI3K), especially the first class of this family of lipid kinases consisting of α, β, γ, and δ subunits, is a promising approach for the treatment of many inflammatory and autoimmune diseases, such as systemic lupus erythematosus or multiple sclerosis [1][2][3]. Because PI3K is involved in many cellular processes, including proliferation, growth, migration, metabolism regulation, and embryogenesis (characterized by high expression of this protein in different human cells), it is considered an excellent therapeutic target [4,5].
Fluorine, which is slightly larger than hydrogen and highly electronegative, constitutes a remarkable role in medicinal chemistry [6][7][8][9][10][11]. Currently, fluorination is a standard strategy to improve the bioavailability of designed drugs and optimize their biological activity [12]. In the last decade, nearly 30% of the drugs approved by the US Food and Drug Administration (FDA) contained fluorine, and fluorinated pharmaceuticals accounted for over 50% of the most profitable drugs worldwide [6,13].
Fluorine, which is slightly larger than hydrogen and highly electronegative, constitutes a remarkable role in medicinal chemistry [6][7][8][9][10][11]. Currently, fluorination is a standard strategy to improve the bioavailability of designed drugs and optimize their biological activity [12]. In the last decade, nearly 30% of the drugs approved by the US Food and Drug Administration (FDA) contained fluorine, and fluorinated pharmaceuticals accounted for over 50% of the most profitable drugs worldwide [6,13].
Molecular modeling methods, such as molecular docking, molecular dynamics (MD), or free-energy perturbation, are widely applied during the rational designing of new compounds to evaluate the formation of ligand-receptor complexes [14]. However, the standard (rigid) docking method has two main limitations: (i) the conformation space is reduced due to the limitations imposed on the system (rigid receptor) and (ii) conventional scoring functions do not consider inductive or resonance effects (which is extremely important for fluorine derivatives). Nevertheless, this method can be successfully used to quickly score poses and find promising hits from a large library of compounds [15,16]. Recently, we proposed a method that can overcome some limitations of the standard molecular docking approach [7,9]. We demonstrated that a workflow comprising a combination of more sophisticated methods, such as induced-fit docking (IFD), quantum polarized ligand docking (QPLD), and binding-free energy calculations based on the Generalized Born Surface Area (GBSA), is more accurate in the prediction of the ligand-receptor complex and its energy than the standard docking procedure, but it is more computationally expensive [17].
In this work, we used the previously described in silico workflow [7,9] to design and determine the biological activity of new fluorinated analogs of CPL302415 ( Figure 1), a PI3Kδ inhibitor. The usefulness of this workflow was validated by correlating the biological activity (IC50) with energy changes calculated in the GBSA model.

Results and Discussion
Molecular docking of CPL302415 derivatives was performed on the previously described PI3Kδ crystal structure (PDB ID: 2WXL). The binding mode observed for CPL302415 (2) was consistent with the previously reported ones [18,19] (Figure 2). The nitrogen of the (difluoromethyl)-1H-benzimidazole fragment interacted with the positively charged Lys779, while the morpholine ring at position 7 (which is required for interaction with PI3Kδ at its catalytic site) formed a hydrogen bond with the main chain of Val828. The tert-butyl piperazine moiety bound to the Trp760 shelf through C-H⋯ π/cation⋯π interaction ( Figure 2).

Results and Discussion
Molecular docking of CPL302415 derivatives was performed on the previously described PI3Kδ crystal structure (PDB ID: 2WXL). The binding mode observed for CPL302415 (2) was consistent with the previously reported ones [18,19] (Figure 2). The nitrogen of the (difluoromethyl)-1H-benzimidazole fragment interacted with the positively charged Lys779, while the morpholine ring at position 7 (which is required for interaction with PI3Kδ at its catalytic site) formed a hydrogen bond with the main chain of Val828. The tert-butyl piperazine moiety bound to the Trp760 shelf through C-H· · · π/cation· · · π interaction ( Figure 2).
According to molecular docking, the static nature of biomolecules is the main source of their limitations, as it does not consider the dynamic nature of biological structures [20]. The analysis of the structure-activity relationship (SAR) of our library revealed that the change from the difluoro group (2) to the methyl group (1) caused a decrease in the activity of the compound, but the introduction of the trifluoromethyl group (3) almost inactivated the compound (Table 1). Interestingly, the docking scores of the above-mentioned derivatives indicated that compound 1 should have the highest IC 50 value (Table 1) in the series-even if we took into account the average docking score for the top three poses, the tendency was almost identical. We substituted position 3 with chlorine (4) or bromine (5) and the resulting compounds were the most potent based on docking scores (Table 1), but they had worse IC 50 values than compound 2.   (2) in the catalytic center of PI3Kδ obtained by molecular docking [19].
According to molecular docking, the static nature of biomolecules is the main source of their limitations, as it does not consider the dynamic nature of biological structures [20]. The analysis of the structure-activity relationship (SAR) of our library revealed that the change from the difluoro group (2) to the methyl group (1) caused a decrease in the activity of the compound, but the introduction of the trifluoromethyl group (3) almost inactivated the compound (Table 1). Interestingly, the docking scores of the abovementioned derivatives indicated that compound 1 should have the highest IC50 value (Table 1) in the series-even if we took into account the average docking score for the top three poses, the tendency was almost identical. We substituted position 3 with chlorine (4) or bromine (5) and the resulting compounds were the most potent based on docking scores (Table 1), but they had worse IC50 values than compound 2.
To assess the accuracy of the presented method, correlation coefficients between the IC50 values and the obtained docking scores were calculated (for the first pose and the average of the first three poses, respectively). The tests showed that there was no correlation (0.53; p > 0.05, and 0.51; p > 0.05).   (2) in the catalytic center of PI3Kδ obtained by molecular docking [19].   (2) in the catalytic center of PI3Kδ obtained by molecular docking [19].
According to molecular docking, the static nature of biomolecules is the main source of their limitations, as it does not consider the dynamic nature of biological structures [20]. The analysis of the structure-activity relationship (SAR) of our library revealed that the change from the difluoro group (2) to the methyl group (1) caused a decrease in the activity of the compound, but the introduction of the trifluoromethyl group (3) almost inactivated the compound (Table 1). Interestingly, the docking scores of the abovementioned derivatives indicated that compound 1 should have the highest IC50 value (Table 1) in the series-even if we took into account the average docking score for the top three poses, the tendency was almost identical. We substituted position 3 with chlorine (4) or bromine (5) and the resulting compounds were the most potent based on docking scores (Table 1), but they had worse IC50 values than compound 2.
To assess the accuracy of the presented method, correlation coefficients between the IC50 values and the obtained docking scores were calculated (for the first pose and the average of the first three poses, respectively). The tests showed that there was no correlation (0.53; p > 0.05, and 0.51; p > 0.05). To assess the accuracy of the presented method, correlation coefficients between the IC 50 values and the obtained docking scores were calculated (for the first pose and the average of the first three poses, respectively). The tests showed that there was no correlation (0.53; p > 0.05, and 0.51; p > 0.05).
Based on the binding mode of CPL302415 (2), we decided to substitute tert-butyl moiety with a benzyl fragment to increase the number of π-π or hydrophobic interactions with Thr750. Therefore, a series of fluorine derivatives were synthesized using the 7-(morpholin-4-yl)pyrazolo[1,5-a]pyrimidin-5-yl]-2-(difluoromethyl)-1H-benzimidazole core (6-9) as well as with an additional carbonyl group in the 1-[2-(4-benzylpiperazine-1-carbonyl) fragment (substituent in position 2 of the pyrazolo[1,5-a]pirymidine core) (10-13) ( Table 2). We docked both nonfluorinated (in the aromatic ring) cores with benzyl fragment compounds to the 2WXL crystal structure, and the observed binding modes were similar to those of CPL302415. However, the hydrogen bonds formed by the morpholine fragment in compound 10 with Val828 were characterized by the worst geometric parameters in comparison to compound 6 ( Figure 3), which could be associated with the lower activity of compound 10 ( Table 2).   Based on the binding mode of CPL302415 (2), we decided to substitute tert-butyl moiety with a benzyl fragment to increase the number of π-π or hydrophobic interactions with Thr750. Therefore, a series of fluorine derivatives were synthesized using the 7-(morpholin-4-yl)pyrazolo[1,5-a]pyrimidin-5-yl]-2-(difluoromethyl)-1H-benzimidazole core (6-9) as well as with an additional carbonyl group in the 1-[2-(4-benzylpiperazine-1carbonyl) fragment (substituent in position 2 of the pyrazolo[1,5-a]pirymidine core) (10-13) ( Table 2). We docked both nonfluorinated (in the aromatic ring) cores with benzyl fragment compounds to the 2WXL crystal structure, and the observed binding modes were similar to those of CPL302415. However, the hydrogen bonds formed by the morpholine fragment in compound 10 with Val828 were characterized by the worst geometric parameters in comparison to compound 6 ( Figure 3), which could be associated with the lower activity of compound 10 ( Table 2).    The carbonyl group affected the orientation of the benzylpiperazine fragment, causing a change in the cation· · · π interaction with Trp760 ( Figure 3). The analysis of the IC 50 values did not allow us to draw any conclusions because in the first series of derivatives (6)(7)(8)(9), the meta-fluoro derivative 8 had the highest IC 50 value, while in the second series (10)(11)(12)(13), the meta-fluoro derivative 12 was the best (Table 2), which suggested that carbonyl oxygen had an impact on binding. Interestingly, we found a poor correlation between any docking score (nor the first best and the best three poses) and biological activity (IC 50 ) (correlation coefficients were −0.02 and 0.03, respectively).
Due to the low correlation of the docking scores obtained in the standard molecular docking approach (Tables 1 and 2) with biological activity, a previously described and validated computational workflow [7,9] was used to compensate for the limitations of conventional scoring functions and increase the accuracy of the pose prediction, especially for fluorinated derivatives ( Figure 4). Because fluorine is the most electronegative element, its substitution leads to significant changes in the distribution of electron density and as a consequence, resonance and inductive effects. The standard molecular docking approach does not consider these effects, while the QPLD method uses atomic charges of a ligand calculated using the quantum mechanical (QM)/molecular mechanical (MM) approach in the protein environment for docking ( Figure 4). To minimize the uncertainty of predicted poses, the binding-free energy was calculated for three poses obtained at the QPLD stage with the smallest root square mean deviation (RMSD) to the core. The energy change (∆∆G) after fluorination was estimated in comparison to the nonfluorinated compound ( Figure 4). Additionally, by using the MD of a nonfluorinated compound, our approach enabled us to relax the binding pocket and fit it into the particular core. Due to its small size [6,21], we assumed that the introduction of a fluorine atom would not result in large conformational changes. The carbonyl group affected the orientation of the benzylpiperazine fragment, causing a change in the cation⋯π interaction with Trp760 ( Figure 3). The analysis of the IC50 values did not allow us to draw any conclusions because in the first series of derivatives (6-9), the meta-fluoro derivative 8 had the highest IC50 value, while in the second series (10)(11)(12)(13), the meta-fluoro derivative 12 was the best (Table 2), which suggested that carbonyl oxygen had an impact on binding. Interestingly, we found a poor correlation between any docking score (nor the first best and the best three poses) and biological activity (IC50) (correlation coefficients were −0.02 and 0.03, respectively).
Due to the low correlation of the docking scores obtained in the standard molecular docking approach (Tables 1 and 2) with biological activity, a previously described and validated computational workflow [7,9] was used to compensate for the limitations of conventional scoring functions and increase the accuracy of the pose prediction, especially for fluorinated derivatives ( Figure 4). Because fluorine is the most electronegative element, its substitution leads to significant changes in the distribution of electron density and as a consequence, resonance and inductive effects. The standard molecular docking approach does not consider these effects, while the QPLD method uses atomic charges of a ligand calculated using the quantum mechanical (QM)/molecular mechanical (MM) approach in the protein environment for docking ( Figure 4). To minimize the uncertainty of predicted poses, the binding-free energy was calculated for three poses obtained at the QPLD stage with the smallest root square mean deviation (RMSD) to the core. The energy change (ΔΔG) after fluorination was estimated in comparison to the nonfluorinated compound ( Figure 4). Additionally, by using the MD of a nonfluorinated compound, our approach enabled us to relax the binding pocket and fit it into the particular core. Due to its small size [6,21], we assumed that the introduction of a fluorine atom would not result in large conformational changes. Computational workflow used to predict the most potent fluorine derivative. The workflow began with the IFD of a nonfluorinated compound (core) (1) followed by 100 ns-long MD simulations (2), which were then clustered based on the RMSD matrix of the protein backbone (3). All derivatives were docked to the three most frequently observed conformations of the protein using the QPLD algorithm (4). Using the MM-GBSA approach, the binding energy (ΔG) was calculated for the three conformations of the ligand with the smallest RMSD of the core to nonfluorinated compounds (5). Finally, the difference in the interaction energy between the most active compound and subsequent isomers (ΔΔG) was calculated (6).
Analysis of the MD trajectories showed that the previously described binding mode [19] and the formed interactions were highly stable ( Figure 5). We found that Lys779  (2), which were then clustered based on the RMSD matrix of the protein backbone (3). All derivatives were docked to the three most frequently observed conformations of the protein using the QPLD algorithm (4). Using the MM-GBSA approach, the binding energy (∆G) was calculated for the three conformations of the ligand with the smallest RMSD of the core to nonfluorinated compounds (5). Finally, the difference in the interaction energy between the most active compound and subsequent isomers (∆∆G) was calculated (6).
Analysis of the MD trajectories showed that the previously described binding mode [19] and the formed interactions were highly stable ( Figure 5). We found that Lys779 involved in a hydrogen bond with benzimidazole nitrogen and tryptophane, which interacted with the positively charged nitrogen of the piperazine group, had high stability and a permanent position. Additionally, Tyr813 was engaged in a weak hydrogen bond with the CH donor ( Figure 5). The substitution of the tert-butyl fragment with the benzyl fragment allowed an additional C-H· · · π/hydrophobic interaction with Thr750. involved in a hydrogen bond with benzimidazole nitrogen and tryptophane, which interacted with the positively charged nitrogen of the piperazine group, had high stability and a permanent position. Additionally, Tyr813 was engaged in a weak hydrogen bond with the CH donor ( Figure 5). The substitution of the tert-butyl fragment with the benzyl fragment allowed an additional C-H⋯π/hydrophobic interaction with Thr750. The comparison of the first top pose for each compound in the series obtained by applying the standard approach showed that the 7-(morpholin-4-yl)pyrazolo[1,5a]pyrimidin-5-yl]-2-(difluoromethyl)-1H-benzimidazole core had less mobility, whereas the tert-butyl 1-{2-[(4-tert-butylpiperazin-1-yl) and 1-[2-(4-benzylpiperazine-1-carbonyl) fragment was directed toward the solvent space, and thus had higher flexibility ( Figure  6I). Since fluorine is a bioisostere of hydrogen [6,8], the substitution of this element should not lead to drastic changes in the binding mode. The comparison of the observed binding modes showed that there was no coherence ( Figure 6I) and that small molecular changes highly affected the position of adjacent fragments of molecules.
The results obtained using different methods of atomic charge assignment are clearly different ( Figure 6I,II). Compared to standard docking, in which the tert-butyl piperazine or benzylpiperazine fragment resulted in high flexibility and mobility ( Figure 6I), the binding poses obtained by the proposed workflow (where the atomic charges are QMderived) showed lower RMSD (less than 0.5 Å; Figure 6II). In the case of the QPLD-docked poses, the introduced fluorine atom did not induce any drastic conformational changes ( Figure 6II). However, due to its high electronegativity, this element may affect the tuning of atomic charges, pKa, and the basicity or acidity (electron density distribution in general) of neighboring functional groups, rather than intermolecular interactions with the biological target [6]. The comparison of the first top pose for each compound in the series obtained by applying the standard approach showed that the 7-(morpholin-4-yl)pyrazolo[1,5-a]pyrimidin-5-yl]-2-(difluoromethyl)-1H-benzimidazole core had less mobility, whereas the tert-butyl 1-{2-[(4-tert-butylpiperazin-1-yl) and 1-[2-(4-benzylpiperazine-1-carbonyl) fragment was directed toward the solvent space, and thus had higher flexibility ( Figure 6I). Since fluorine is a bioisostere of hydrogen [6,8], the substitution of this element should not lead to drastic changes in the binding mode. The comparison of the observed binding modes showed that there was no coherence ( Figure 6I) and that small molecular changes highly affected the position of adjacent fragments of molecules.
The results obtained using different methods of atomic charge assignment are clearly different ( Figure 6I,II). Compared to standard docking, in which the tert-butyl piperazine or benzylpiperazine fragment resulted in high flexibility and mobility ( Figure 6I), the binding poses obtained by the proposed workflow (where the atomic charges are QM-derived) showed lower RMSD (less than 0.5 Å; Figure 6II). In the case of the QPLDdocked poses, the introduced fluorine atom did not induce any drastic conformational changes ( Figure 6II). However, due to its high electronegativity, this element may affect the tuning of atomic charges, pK a , and the basicity or acidity (electron density distribution in general) of neighboring functional groups, rather than intermolecular interactions with the biological target [6]. Using the MM-GBSA approach, the interaction energy of the obtained ligandreceptor complex was calculated as the average of the top three ligand poses (∆ ) (Tabs. 3 and 4). The correlation coefficient of the biological activity (IC50) and the MM-GBSA scores of the tert-butyl piperazine derivatives was significantly higher with the use of our approach compared to standard docking scores (0.95; p < 0.05 and ~0.5, respectively) ( Table 3). The results obtained for the tert-butyl piperazine derivatives showed that the introduction of a difluoro group to the 2-methyl-1H-benzimidazole fragment increased the stabilization energy of the complex (ΔΔG = −7.6), whereas a trifluoro group decreased the stability of complexes, which was probably due to steric hindrances and inductive effects (Table 3). Moreover, as shown by Kurczab et al. [22], the MM-GBSA approach can be successfully used for heavier halogens. Therefore, we extended our library with derivatives containing bromine or chlorine at position 3. These derivatives were less active, but the proposed algorithm accurately predicted the energy loss for the difluoromethyl-1H-benzimidazole derivative (Table 3). For 4-benzylpiperazine and 4-benzylpiperazine-1-carbonyl fluorine derivatives, the correlation coefficient was almost equal to 1 (0.9; p < 0.05 and 0.95; p < 0.05, respectively), which implies that although there was no correlation with the docking scores, we obtained a perfect correlation here (Table 4).  (Tables 3 and 4). The correlation coefficient of the biological activity (IC 50 ) and the MM-GBSA scores of the tert-butyl piperazine derivatives was significantly higher with the use of our approach compared to standard docking scores (0.95; p < 0.05 and~0.5, respectively) ( Table 3). The results obtained for the tert-butyl piperazine derivatives showed that the introduction of a difluoro group to the 2-methyl-1H-benzimidazole fragment increased the stabilization energy of the complex (∆∆G = −7.6), whereas a trifluoro group decreased the stability of complexes, which was probably due to steric hindrances and inductive effects (Table 3). Moreover, as shown by Kurczab et al. [22], the MM-GBSA approach can be successfully used for heavier halogens. Therefore, we extended our library with derivatives containing bromine or chlorine at position 3. These derivatives were less active, but the proposed algorithm accurately predicted the energy loss for the difluoromethyl-1H-benzimidazole derivative (Table 3).  For 4-benzylpiperazine and 4-benzylpiperazine-1-carbonyl fluorine derivatives, the correlation coefficient was almost equal to 1 (0.9; p < 0.05 and 0.95; p < 0.05, respectively), which implies that although there was no correlation with the docking scores, we obtained a perfect correlation here (Table 4).
In addition, to demonstrate the superiority of the proposed workflow, the QPLD docking protocol was performed on every stage of kinase flexibility: (i) on rigid crystalized PI3Kδ structure (PDB ID: 2WXL) and (ii) on the induced-fit docked poses chosen for the MD stage. The molecular dynamics were the most time-consuming step; therefore, if, thanks to QM-derived atomic charges (which could help obtain a more accurate ligandreceptor interaction energy), a higher correlation of IC 50 values with ∆G values in the previous stages could be achieved, the computational time could be saved compared to the presented workflow. In the same way as in the presented workflow, using the MM-GBSA approach, the interaction energy of the obtained ligand-receptor complex was calculated as the average of the top three ligand poses ( − ∆G) (SI Tables S1-S4). For the rigid conformation of the catalytic pocket, the correlation coefficients of the biological activity (IC 50 ) and the MM-GBSA score of tert-butyl piperazine derivatives were slightly lower (−0.34; p > 0.05) compared to the standard docking scores (~0.5) and proposed workflow (0.95) (SI Table S1). The correlation coefficients for 4-benzylpiperazine (−0.43; p > 0.05) and 4-benzylpiperazine-1-carbonyl fluorine derivatives (−0.78; p > 0.05) were better than this obtained with the standard (rigid) approach (~0.0) but lower than those obtained with our workflow (~0.93). Next, the QPLD approach was used on grids obtained in the IFD protocol for compounds 1, 6, and 10. The correlation coefficient for tert-butyl piperazine derivatives was similar (−0.35; p > 0.05) to that obtained with a rigid catalytic pocket (−0.34), whereas the correlation was better for 4-benzylpiperazine (0.60; p > 0.05) and 4-benzylpiperazine-1-carbonyl fluorine derivatives (−0.78; p > 0.05) compared to those obtained with the non-flexible structure of PI3Kδ (−0.43; −0.78, respectively). Nevertheless, the correlation coefficient obtained with the described workflow was higher (~0.93), which confirms the importance of all stages in the proposed workflow.
It is worth stressing that the correlation of IC 50 values with ∆G values increases with greater flexibility of the catalytic pocket and its fit to the molecular core. These results suggest that properly prepared catalytic pockets for chemical cores in the IFD and MD stages, combined with QM-derived atomic charges, can be effectively used for prediction as well as to improve the discrimination of active compounds from inactive ones.

Chemistry
The compounds 1-14 discussed in this work were synthesized following the general procedures presented in our previous papers [19,23]. The synthesis pathway has multiple steps, of which the last one is described here. The final compounds were synthesized via reactions such as reductive amination, amidation, or coupling reactions including the Buchwald-Hartwig reaction, and their yields varied depending on the structure. The Buchwald-Hartwig reaction, amidation reaction, and reductive amination are especially described in this work.

General Information
The reagents (at least 95% purity) were purchased from ABCR (Dallas, TX, USA), Acros (Geel, Belgium), Alfa Aesar (Haverhill, MA, USA), Combi-Blocks (San Diego, CA, USA), Fluorochem (Hadfield, UK), (Buchs, Switzerland), Merck (Darmstadt, Germany), and Sigma Aldrich (Saint Louis, MI, USA) and were used without additional purification. Solvents were purified according to the standard procedures if required. Air-or moisture-sensitive reactions were carried out under an argon atmosphere. The progress of all reactions was routinely monitored by thin-layer chromatography (TLC). TLC was performed on silica gel-coated plates (Kieselgel F254), which were visualized using a UV light. Flash chromatography was performed on Merck silica gel 60 (230-400 mesh ASTM). 1 H NMR spectra were acquired using JOEL JNMR-ECZS 400 and 600 MHz spectrometers ( 1 H observed at 400, and 600 MHz, respectively). 13 C NMR spectra were recorded at 101 and 151 MHz, respectively. Due to the poor solubility of some final compounds, the usual characterization using 13 C NMR was omitted. Chemical shifts for 1 H and 13 C NMR spectra were reported in d (ppm) using tetramethylsilane as an internal standard or based on the residual undeuterated solvent signal (2.50 ppm for DMSO-d 6 and 7.26 ppm for CDCl 3 ). The abbreviations for multiplets of 1 H signals are as follows: s (singlet), d (doublet), t (triplet), m (multiplet), dd (doublet of doublets), dt (doublet of triplet), and q (quartet). Coupling constants (J) are expressed in Hertz. A JEOL Royal HFX probe head was used for recording the 13 C NMR spectrum as it allows measurements to be taken with the simultaneous decoupling of both 1 H and 19 F nuclei [24]. Atmospheric pressure ionization and electrospray ionization mass spectra were obtained using an Agilent 6130 LC/MSD spectrometer or Agilent 1290 UHPLC system coupled with an Agilent QTOF 6545 mass spectrometer. All spectra of the final compounds are shown in the Supplementary Materials.

Synthesis
Compounds 1 and 2 were synthesized according to the procedures described in our previous publications [19,23].

General Procedure for the Reductive Amination Reaction
To the solution of the corresponding aldehyde (1.0 eq) in a dry DCM (10 mL/1 g aldehyde), an amine derivative (1.2 eq) was added, and the reaction mixture was stirred for 1 h at room temperature. Then, sodium triacetoxyborohydride (1.5 eq) was added and the mixture was stirred for a further 15 h at room temperature. Next, water was added to the reaction mixture and the phases were separated. The aqueous phase was extracted three times with DCM, while the combined organic phases were dried over anhydrous sodium sulfate, filtered, and concentrated. The resulting residue was purified by flash chromatography.

Computational Workflow Used to Predict the Most Potent Fluorine Derivative
We used a previously described computational workflow involving IFD, MD simulations, and QPLD combined with energy calculations (applying the MM-GBSA method). The crystal structure of PI3Kδ protein (PDB ID: 2WXL) [25] that was successfully used in our previous study to support the SAR analysis [19] was retrieved from Protein Data Bank [26,27].

IFD
The three-dimensional structures of the ligands were prepared using LigPrep v3.6 [28], and appropriate ionization states at pH 7.0 ± 0.5 were assigned using Epik v3.4 [29]. The Protein Preparation Wizard tool [28] was used to assign bond orders and appropriate amino acid ionization states and to check for steric clashes for the PI3Kδ crystal structure. The receptor grid was generated (OPLS3 force field [30]) by centering the grid box with a size of 8 Å on crystalized ligands (ZSTK474). Automated flexible docking [31,32] of nonfluorinated compounds was performed using Glide v6.9 [33][34][35] at the SP level.

MD
MD simulations (100 ns long) were performed using Schrödinger Desmond software [36,37]. Each ligand-receptor complex selected based on the IFD analysis was joined with the POPC (309.5 K) membrane bilayer. The system was solvated by water molecules described by the TIP4P potential [38] in orthorhombic box with a distance of 10Å from the complex, using the OPLS3e force field [30] for all atoms. NaCl (0.15 M) was added to mimic the ionic strength inside the cell. The simulations were carried out using the NPAT protocol at a temperature of 309.5 K and a pressure of 1013.25 hPa.

QPLD
The grids for the receptors were generated (OPLS3e force field) by centering the grid box (size 8 Å) on a ligand. Docking of all fluorinated compounds was performed by the QPLD [16] procedure involving the QM-derived ligand atomic charges in the protein environment at the BLYP/cc-pVDZ level [39,40]. Five poses were obtained for each ligand.

Binding-Free Energy Calculations
Using GBSA, the binding-free energy was calculated based on the ligand-receptor complexes generated by the QPLD procedure. The ligand poses were minimized using the local optimization feature in Prime. The distance of the flexible residue from a ligand pose was set to 6 Å. Ligand charges obtained in the QPLD stage were used. The energies of complexes were calculated with the OPLS3e force field and the GBSA continuum solvent model. To assess the influence of a given substituent on binding, ∆∆G was calculated as the difference between the binding-free energy (∆G) of the nonfluorinated compound and its fluorinated analogs.

Conclusions
We evaluated a workflow involving IFD, MD, and QM/MM docking (QPLD) with the MM-GBSA calculation to score halogenated derivatives of CPL302415, a clinical PI3Kδ selective inhibitor, and compared the results with those obtained by the Glide scoring function. Additionally, we synthesized a series of novel fluorinated compounds to estimate the accuracy of the pose prediction in the molecular docking procedure.
We found that a properly prepared catalytic (binding) pocket for chemical cores in the IFD and MD stages, combined with QM-derived atomic charges, can be effectively used for prediction, as well as to improve the discrimination of active compounds from inactive ones. Moreover, the standard approach (rigid docking) seems to be less effective at scoring halogenated derivatives due to the fixed atomic charges and the fact that it does not consider the response and indictive effects caused by fluorine. Our results suggest that the proposed computational workflow may be a valuable tool for the rational design of new halogenated drugs.