Wolfram Cloud Document

Interleukin proteins present vast opportunities for treating autoimmune and inflammatory diseases through proper regulation. While antibody treatments exist for some interleukins, their expensive manufacturing costs prevent their widespread use. Numerous drug screenings for small-molecule inhibitors have been performed for a variety of interleukins, but in-depth quantum chemistry computations are yet to be applied to these candidates to understand their binding properties. In this exploration, quantum and molecular mechanics techniques for energy optimization are applied to two potential interleukin inhibitors and their respective interleukin proteins. A method to compute and model small molecule binding with proteins in Wolfram is created along with a way to visualize residues and corresponding atom indices in a Molecule object. From the binding of the small-molecule lanraplenib to interleukin-8, lanraplenib was found to bind to the same binding site of interleukin-8 to its receptor, potentially acting as a interleukin-8 inhibitor.

Introduction

Interleukin proteins are a group of cytokines responsible for cell signaling and communication across the immune system and body cells. Interleukins encompass a vast range of essential roles for the immune system, from activation, differentiation, proliferation, migration, and adhesion. This group of small proteins are secreted by immune cells and are primarily synthesized by CD4 helper T-lymphocytes, macrophages, monocytes, and epithelial cells, due to their roles as initial messengers for the immune system. Once secreted, interleukins bind to a variety of specific receptors on cell surfaces and regulate the immune response; their varied functions can elicit up-regulatory and down-regulatory responses in cells, leading to both pro-inflammatory and anti-inflammatory effects (Justiz Vaillant & Qurie, 2025). However, interleukins’ large impact on the immune system means abnormal production and regulation of certain interleukins are linked to a variety of autoimmune diseases, such as lupus, rheumatoid arthritis, psoriasis. Certain interleukins’ pro-inflammatory functions play a role in inflammatory diseases, such as chronic obstructive pulmonary disease (COPD) and asthma. In addition, interleukins can be utilized in immunotherapy and cancer treatments, where either the promotion or inhibition of their production can be promising treatments (Akdis et al., 2016).

Identifying a small-molecule inhibitor to certain interleukins, like interleukin-4 (IL-4), interleukin-6 (IL-6), and interleukin-8 (IL-8) is a promising alternative to current inhibitory treatments, which require the expensive and lengthy synthesis of specialized antibodies. Compared to large antibodies, small molecules can more easily penetrate tissues and have cheaper production costs. However, identifying an effective ligand that prevents the protein-protein interaction between interleukins and their receptors is a challenging endeavor, as the necessary chemotypes needed for proper inhibition have not been thoroughly identified, and only a handful of interleukins have identified small-molecule inhibitors (Quinnell et al., 2020).

Drug discovery for a small-molecule inhibitor for interleukins are primarily carried out in silico through pharmacophore model generation, docking, and molecular dynamics simulations to find drug candidates. These computational algorithms rely on molecular mechanics (MM) for their calculations on molecule binding and structure stability. Molecular mechanics uses principles of classical mechanics to approximate molecule interaction. In quantum mechanics (QM), the Schrödinger equation defines molecule interaction by defining the movement atomic particles. Solving the Schrödinger equation requires heavy computations, and there are varying approximations and methods to solve this equation. Molecular mechanics does not solve this equation, and does not calculate the fine details of molecule/atomic interaction, leading to significantly faster computations. While quantum mechanics is suitable for small molecules, proteins like interleukins have significantly more atoms to account for.

A quantum mechanics/molecular mechanics (QM/MM) approach to modeling protein-small molecule interaction combines both accuracy and speed to analyze the specific chemical interactions between a small molecule and larger molecule, like proteins. A molecular mechanics optimization with a molecular force field uses a generalized set of parameters to approximate potential energy of the current structure, and systematically changes coordinates and molecular properties until the potential energy is minimized. A force field’s generalizations of molecule interaction results in fast computation times, which is important with protein optimization. Combined with a quantum mechanics optimization for the small molecule interaction with a specific residue or group of residues, more accurate information on the specific chemical interactions between the small molecule and the binding site can be collected. A quantum mechanics/molecular mechanics approach to modeling interleukin protein-small molecule interaction is presented, with both amino acid-small molecule quantum mechanics calculations and protein-small molecule molecular mechanics calculations.

Modeling Proteins

In a 2021 docking study, Prathap et al. identified binding of stachydrine, a small-molecule that naturally occurs in citrus, to interleukin-6. It was found that stachydrine forms a hydrogen bond to glutamic acid (Prathap et al., 2021). A replication of these results are simulated in Wolfram, and quantum chemistry computations were run to analyze the glutamic acid-stachydrine interaction further. Before the molecular/quantum mechanics calculations can be performed, the molecules must be prepared first.

Setting Up Molecules

Build stachydrine from its IUPAC name:

In[]:=

molStachydrine=Molecule["(2S)-1,1-dimethylpyrrolidin-1-ium-2-carboxylate"]

Out[]=

Molecule

Formula:

Atoms: 23 Bonds: 23



Experimental data collected on Interleukin-6 can be downloaded and imported from the Protein Database (PDB), with the PDB ID “1IL6” (Xu et al., 1997).

Pulling in the experimental data using ExternalIdentifier[] with the PDB Structure ID:

In[]:=

interleukin6BioMolecule=BioMolecule[ExternalIdentifier["PDBStructureID","1IL6"]]

Out[]=

BioMolecule

Type: Peptide

Chains: 1



To work with the molecules, interleukin-6 must be converted into a Molecule object. Since BioMolecue does not explicitly model hydrogens, the hydrogens must be made explicit when converted.

Convert the BioMolecule into a Molecule object and make hydrogens explicit:

In[]:=

interleukin6Mol=Quiet@MoleculeModify[Molecule[interleukin6BioMolecule],"MakeHydrogensExplicit"]

Out[]=

Molecule

Formula:

832

1372

229

261

Atoms: 2703 Bonds: 2720



Cleaning Protein Data

Proteins imported from the protein database is derived from experimental data using techniques such as X-ray crystallography, solution and solid-state nuclear magnetic resonance (NMR) spectroscopy, and cryogenic electron microscopy (cryo-EM). The vast majority of these techniques require an aqueous protein sample, and the polar nature of water and proteins means water molecules typically bind to the protein, both on its surface and in its internal cavities.

Extra molecules not belonging to the protein, like the observed water molecules are retained in this experimental data. Some proteins from the database have these extra molecules recorded as a separate unit that can be cut out from the BioMolecule[] object. If the data has been separated, BioMolecule will automatically identify the different chain-types which can then be cleaved out. The protein model of interleukin-6 used (1IL6) does not contain water molecules, but the majority of protein data will contain water molecules. Sometimes, non-polymer ligands, like sulfates, will also be contained in the protein data.

The interleukin-6 protein data used previously was obtained through solution NMR. The X-ray crystallography protein data of interleukin-6 (1ALU) contains water molecules and additional small molecules that need to be cleaned before using in computations (Somers et al., 1997).

Import interleukin-6 data using the PDB Structure ID:

In[]:=

interleukin6XRay=BioMolecule[ExternalIdentifier["PDBStructureID","1ALU"]]

Out[]=

BioMolecule

Type: Peptide

Chains: 7



Grab the chain-type information of IL-6:

In[]:=

interleukin6XRay["ChainTypes"]

Out[]=

APeptide,BNonPolymer,CNonPolymer,DNonPolymer,ENonPolymer,FNonPolymer,GWater

Notice that there are 7 chains in the BioMolecule object, with chain A being the actual peptide chain and Chain G being the group of water molecules. Chains B-G can be removed from the BioMolecule, resulting in just the isolated peptide.

BioMoleculePlot3D[] does not represent water molecules in its plotting, but it does include any non-polymers present in the data.

Plot the uncleaned X-ray crystallography protein:

In[]:=

BioMoleculePlot3D[interleukin6XRay,ImageSize->Large]//Rasterize

The molecular structure of the protein will include any water molecules floating around or bound to the protein from the data.

Visualize the floating water molecules in orange:

In[]:=

MoleculePlot3DMolecule[interleukin6XRay],StyleMoleculePattern

Atoms: 3

Bonds: 2

,Orange//Rasterize

Out[]=

To remove both non-polymers and waters from the protein, the peptide chain should be isolated.

To isolate chain A, form a new BioMolecule object by grabbing the peptide object

In[]:=

dehydratedInterleukin6XRay=BioMolecule[interleukin6BioMolecule-><|"Chains"->{"A"}|>]

Out[]=

BioMolecule

Type: Peptide

Chains: 1



Visualize the cleaned protein in the same way as the solution NMR IL-6:

In[]:=

BioMoleculePlot3D[dehydratedInterleukin6XRay,ColorRules->"Residues",ImageSize->Large]//Rasterize

The BioMolecule plot is identical to the original version of the protein as BioMoleculePlot3D[] only registers the peptide compounds in the plotting.

Visualize the Molecule plot of the protein:

In[]:=

MoleculePlot3D[Molecule[dehydratedInterleukin6XRay],PlotTheme->"Tubes",ImageSize->Large]//Rasterize

Out[]=

Energy Calculations

A Merck molecular force field (MMFF) energy calculation runs a single point calculation of a given molecule in its current formation. Treating the IL-6 protein and the stachydrine small molecule as a singular molecule tests its energy in relation to each other. Multiple single point calculations with modifications to the two molecules’ locations and bonds optimizes its geometry, and an optimized geometry is equivalent to a simulated binding site.

Treating the IL-6 protein and the stachydrine as a combined molecule requires the MoleculeComplex paclet, which, as it’s name suggests, forms a molecule complex that combines the two molecules as one object.

Install the MoleculeComplex paclet, which can be accessed through the paclet repository.

A MoleculeComplex[] object requires all of its input molecules to have 3-dimensional coordinates. The Molecule object of IL-6 has 3-dimensional coordinates defined already by the database import. Stachydrine is only defined as a 2-dimensional structure currently.

Use MakeMolecule3D[] from the MoleculeComplex paclet to calculate stachydrine’s 3-dimensional coordinates

MoleculeComplex[] returns a MoleculeComplex object made up of its input molecules, placing them together.

This MoleculeComplex can be used to run molecular mechanics computations with the Merck molecular force field shown later with interleukin-8.

Amino Acid & Small Molecule Interaction

Focusing on the specific amino acid-small molecule interaction in the protein is a more manageable and customizable computation, which aids in understanding the quantum chemistry functions before running more expensive, large computations with the whole protein.

Specific amino acids can be pulled either using their 3-letter or single letter abbreviations

Pull glutamic acid and convert it into a molecule:

Perform the same complex operation with stachydrine and the amino acid.

MoleculeComplex currently just overlays the stachydrine and amino acid together into the same object without any optimization.

Of course, this molecule configuration can only exist in simulation. The single point energy calculation of the MoleculeComplex’s current configuration gives some indication on the stability of the compound.

MoleculeComplexEnergy[] can calculate the MMFF-based single-point energy of the molecule at its current configuration:

The single point energy calculation of this current configuration is absurdly high, and definitely does not occur naturally. Finding the naturally optimal interaction between stachydrine and glutamic acid requires geometry optimization by calculating electron structures of the proteins.

Energy & Geometry Optimization of Molecule Complexes

Running a geometry optimization will iteratively test a wide variety of configurations to find the lowest possible energy and simulate a more accurate representation of the molecules together.

Use MoleculeComplexOptimizeGeometry[] to run a MMFF-based wave function approximation to find the optimal geometry of the molecule complex:

Checking the new MMFF-based single-point energy after optimization through the same process as earlier:

The negative single-point energy, compared to the earlier’s energy, indicates a much more stabilized structure, and signals the possibility of newly formed intermolecular forces. Hydrogen bonds formed between the stachydrine and glutamic acid are calculated by measuring the bond lengths of this optimized molecule interaction.

Plot the optimized geometry using MoleculeComplexPlot3D[]. Label the atom indices to find out which molecules to test bond lengths through the “AtomLabels” option.

The hydrogen bonds that can be formed between the stachydrine and glutamic acid are formed by the oxygens’ negative partial charge. Atoms of index 18, 8, and 4 of glutamic acid appear to interact with the atoms of index 9, 15, and 12 of the stachydrine, respectively.

Calculate the bond length between atoms through ComplexMoleculeDistance[]. Molecules in molecule complex are indexed by their molecule index and atom index within the molecule.

Calculations for single point energy and optimizing geometry can also be done in a separate quantum chemistry paclet. While MMFF bases itself on a dataset of previously done calculations/optimizations, it does not give the most precise or accurate results compared to other methods to solve the Schrödinger equation.

Hartree-Fock and Density-Functional Theory models can better approximate the electronic structure by solving the Schrödinger equation compared to the MMFF models. Specific basis sets can also be defined, but computation becomes much more expensive with these higher level models. The STO-3G basis set is the most simplified set for the Hartree-Fock calculations, and the 6-311G basis set suffices for the most precise computations that will be done between small molecules. For protein-small molecule interaction, MMFF calculations are the highest level optimizations that can be expected for the high atom count of proteins.

Optimized Energy Computations of Stachydrine and Glutamic Acid

Hartree-Fock calculations with basis sets STO-3G and 6-311G, along with a Density Functional Theory calculation using the 6-311G basis set were computed. Bond length differences with each computation are noted. When computing a more precise model or basis set, using the previously computed optimal configuration saves time and resources.

First, import the quantum chemistry paclet into the notebook.

Convert the Molecule Complex object into a Molecule object to run the quantum chemistry on the molecules.

Establish the model to compute the optimal geometry with using ModelChemistry[]

Run the OptimizeMoleculeGeometry[] function on the molecules using the set model, and ensure that the result is assigned to a separate variable.

hfSTO3GStachydrineGlu=OptimizeMoleculeGeometry[moleculeStachydrineGluMMFF,hartreeFockSTO3G]

Previously computed data can be accessed by pulling its information from the cloud.

Visualize the Hartree-Fock STO-3G optimization of stachydrine and glutamic acid by grabbing the Molecule object from the ElectronicStructureResult:

This Molecule object represents the new optimized coordinates. Compare the Hartree-Fock STO-3G optimization’s bond lengths to the MMFF’s optimization to understand their differences.

Repeat the same process done to compute bond distances for MMFF. MoleculeValue is used instead because the computation results in a Molecule object rather than the MoleculeComplex. Atom indices also vary from the conversion of MoleculeComplex to Molecule.

Atom lengths have decreased in this optimization, but this is not necessarily the case for further computations. The Hartree-Fock 6-311G and Density Functional Theory 6-311G computations have also been performed.

Grabbing the Hartree-Fock 6-311G ElectronicStructureResult from the cloud:

Grabbing the Density Functional Theory 6-311G ElectronicStructureResult

Computing the bond lengths of the Hartree-Fock 6-311G and DFT 6-311G computations

Bond lengths are not a perfect tool for comparison. Each model and basis set considers a different range of properties, which lead to these differences in the bond lengths. However, from the results of the two higher level computations indicate that there may not be a hydrogen bond between the stachydrine’s hydrogen (index 12) and glutamic acid’s oxygen (index 27).

Atom Substitution of Small Molecules

From the hydrogen bonds formed between stachydrine and glutamic acid, modifying stachydrine through atom substitution can yield a stronger binding affinity to glutamic acid. The previous optimizations indicate that the bond distance between the double bonded oxygen on the carboxyl group of glutamic acid and the carbon bonded to the nitrogen is too long to form a hydrogen bond. Substituting the carbon atom on stachydrine for a silicon atom could be more favorable with its decreased electronegativity.

Use MoleculeModify to replace an atom based on their indices. The previously optimized molecule structure is used. (CLARIFY WHY)

Run the MMFF optimization of this substitution to determine potential hydrogen bonding. MoleculeModify[] runs an MMFF computation:

Visualize the optimization of the substitution:

Calculate the bond lengths of the optimized silicon substitution:

While the hydrogens connected to the substituted silicon atom are less tightly bonded, the optimized molecular geometry of the silicon-nitrogen molecules separates the potential hydrogen bond further, instead of closer together.

Protein & Small Molecule Interactions

Modeling protein and small molecule ligand interactions can also be computed, although smaller model approximations for the electron structure need to be used for the protein’s large size. Once potential binding sites have been identified for the ligand on the protein, the interacting amino acids can be isolated and analyzed further with more precise computations, as the rest of the protein does not have a significant impact on the ligand’s interaction.

A smaller interleukin protein can be modeled first for its smaller computation sizes. Interleukin-8 (IL-8), a small cytokine composed of only 72 residues after cleaving, is produced primarily by macrophages to induce chemotaxis (directional migration) to a of neutrophils and other white blood cells. IL-8 plays a major role in cell recruitment and immune system activation, but when abnormally regulated, overproduction of IL-8 is linked to a variety of inflammatory diseases. Chronic obstructive pulmonary disease (COPD), asthma, and psoriasis have all been connected to the overexpression of IL-8 (Matsushima et al., 2022).

In a 2022 in silico study by Tran et al., potential small molecule IL-8 inhibitors were identified through a virtual screening of drug databases, using pharmacophore models (using the structure of the target protein to model the way a ligand would bind and cause a biological effect) and later molecular docking on probable candidates (Tran et al., 2022). The protein-ligand interaction with one of the study’s best candidates, lanraplenib (DB14770), can be modeled in Wolfram.

Forming a MoleculeComplex

The imported IL-8 protein needs to be cleaned in the same process as IL-6 before any modifications are performed on it. In addition, some valence errors in the protein are fixed as well. The imported IL-8 protein (PDB ID: “4XDX”) was obtained through X-ray diffraction (Ostrov et al., 2015).

Grab IL-8 information from the protein database:

Remove water molecules by chain isolation:

Grabbing a formal-charge cleaned molecule derived from the dehydratedInter8 BioMolecule:

Import lanraplenib through it’s IUPAC name and calculate 3-D coordinates:

Just like stachydrine and glutamic acid, IL-8 and lanraplenib can be combined into a MoleculeObject together. However, the protein conversion from BioMolecule to Molecule alters the hydrogen atom processing in the Molecule object. BioMolecule objects automatically make hydrogens implicit, essentially hiding them. When converted into a Molecule object to be used with MoleculeComplex, the hydrogens are still implicit, causing the MoleculeComplex function to assume that the protein is missing those hydrogens. Making the hydrogens explicit as a Molecule object overcomes this issue.

Use MoleculeModify[] to make hydrogens explicit in IL-8:

Make a MoleculeComplex combining the IL-8 protein with lanraplenib:

Visualize the combined MoleculeComplex in one plane

To highlight the lanraplenib in a MoleculeComplex, the “MoleculeHighlighting” option is used.

Highlight lanralenib using “MoleculeHighlighting” and specify that all atoms in the molecule should be highlighted:

When combined into a MoleculeComplex, lanraplenib is placed in an open space that does not interfere with IL-8. However, for a more successful energy optimization, lanraplenib needs to be moved closer to potential binding sites on the protein. However, the Molecule object does not store the BioMolecule residue data, meaning a function needs to be made to label atom indices with their associated residue. Once residues have been labeled by their atom index, the lanraplenib can be moved closer to specific residues for optimization. When lanraplenib has been placed into its potential binding site, the molecular mechanics optimization will adjust molecules both in IL-8 and lanraplenib to find the lowest electronic structure energy. The movement of lanraplenib and IL-8 means that lanraplenib may be moved into a completely different binding site than predicted. However, the molecular mechanics optimization does not test moving lanraplenib to completely different coordinates (like on the other side of the protein), unless its optimization leads to that positioning. The molecular mechanics optimization ends once it has reached a relatively more stable position compared to the other nearby areas.

Residue Labeling and Visualization

The BioMolecule object obtained by the protein database contains the ordered list of residues and a list of atomic symbols grouped by their respective residues. When converting to a Molecule object, the atom indices are retained but the residue groupings are not. In addition, some proteins may have hydrogens already included in the database, meaning that they are given indices like any other atom in the backbone, while other proteins will not include hydrogens in the database or BioMolecule object, meaning that their indexes are placed after all other atoms in the backbone. To approach proteins with integrated hydrogen, residues and their indices and residue atom groupings by atom index are needed to form an association linking the residue and its respective atoms. With proteins that have non-integrated hydrogens, the same information as proteins with integrated hydrogens is needed along with integrating connected hydrogen single bonds with their respective residue groups.

Labelling Proteins with Integrated Hydrogens

Gather the list of residues from the BioMolecule and add a label for the residue’s place in the protein:

Since both the BioMolecule atomic symbols and residue lists are ordered, they can be combined and converted into a list of atom indexes.

Convert from atomic element grouping by residue to atom index grouping by residue:

Create an association between residues labeled with their indices and corresponding atom indices:

Labelling Proteins with Non-Integrated Hydrogens

The hydrogen atoms in non-integrated proteins are not present in the BioMolecule’s atomic symbols, as they are implicit. In addition, they will be placed at the end of the Molecule object’s indexing as the conversion into Molecule is when they are added. Because the hydrogens are missing from the BioMolecule object, the prior functions will ignore the hydrogens bonded to the other atoms. By finding the newly added bonds formed by the conversion into a Molecule object, the hydrogen atoms can be added into the created residue groupings by their index. The index threshold for hydrogen bonds is derived by the number of non-hydrogen atoms in the molecule.

Calculate the threshold for the index of hydrogens:

Identify hydrogen bonds, add them to corresponding residue groups, and remove duplicates:

Combine into an association as earlier, using newly defined functions instead:

Visualizing Residues

Once a protein has been identified as having integrated hydrogens or implicit hydrogens, a 3D plot of the protein with desired residues highlighted can be made. Using DynamicModule[] and Manipulate[] creates an interactive interface for a variety of options on the MoleculePlot3D[]. Atom labels can be toggled for the highlighted residues.

Create a function to plot a given protein and highlight residue associations:

In the 2022 study by Tran et al., lanraplenib was predicted to form hydrogen bonds to primarily arginine-6 (ARG-6) and leucine-5 (LEU-5), along with glutamic acid (GLU-4) (Tran et al., 2022). The current protein data used does not model the first 2 residues as it only models 70 out of 72 total residues coded for in the protein, so the labeling of the residues is 2 values less than the complete protein model. The four potential binding residues can be visualized using the plot model. From the visualization, likely hydrogen bond interactions can be used to move lanraplenib to the molecule complex, and specific atom indices belonging to certain residues can be recorded.

First, define the residue associations of IL-8:

Visualize the protein and residues, selecting GLU-2, LEU-3, and ARG-4:

Adjusting a Small Molecule’s Position

To move the lanraplenib to a predicted binding site on IL-8, a temporary bond can be formed between an atom on the desired binding site and an atom on lanraplenib. The temporary bond can be modified using MoleculeModify[], adjusting bond lengths, bond angles, and dihedral (torsion) angles. Alternatively, to better visualize lanraplenib’s position on the protein, MoleculeManipulate3D[] from the Wolfram Function Repository can be used, which creates an interactive screen to adjust lanraplenib.

Convert the Interleukin-8-lanraplenib complex into a molecule so MoleculeModify[] can operate on the complex:

Most interactions between proteins and small molecules are composed of hydrogen bonds. Determining accessible hydrogen bonds that could be formed between IL-8 and lanraplenib serves as a starting point for the molecular mechanics optimization.

Visualize lanraplenib with atom labels to better understand its structure:

From analyzing the potential binding residues and the structure of lanraplenib, a hydrogen bond between the oxygen on GLU-2’s carboxyl group (index 14) and a hydrogen on lanraplenib’s ether group (index 51) seems plausible. In addition, a hydrogen bond between the oxygen on arginine’s carboxyl group (index 43) and a hydrogen on one of lanraplenib’s rings (index 42) is also possible. To move lanraplenib to nest inside of this protein structure, one bond is formed and modified with at a time, as multiple bonds connecting lanraplenib to IL-8 simultaneously results in segmented movement of lanraplenib.

To form a new bond using MoleculeModify, the index of both atoms need to be known. Since lanraplenib was added as the second molecule in MoleculeComplex, it’s molecule indices are directly after IL-8’s. Since IL-8 has 1165 atoms, 1165 is added to lanraplenib’s atom indices.

Create a temporary hydrogen bond between IL-8 and lanraplenib:

Manipulate the bond properties of the temporary bonds to move lanraplenib:

The adjusted molecule after adjusting bond lengths and angles can be saved as a new molecule object:

Plot the adjusted IL-8-lanraplenib complex, with a temporary hydrogen bond formed between the ether group and glutamic acid:

Remove the temporary hydrogen bond:

To better align lanraplenib to IL-8, another adjustment with a temporary hydrogen bond is done.

Repeat the same process by forming a different bond:

Adjust the bond length and angles:

Save the adjusted molecule:

Plot the adjusted molecule:

Currently, this modification moves the lanraplenib into IL-8 -- this will be fixed during the energy optimization.

Remove the temporary hydrogen bond:

This structure definitely cannot exist naturally as there are some atoms that are very close together. The optimization of this molecule has been performed already.

Notice that lanraplenib (in purple) has shifted positions. Instead of originally being placed on the side, it has now docked into a valley-shaped portion of IL-8. This structure has the lowest energies out of all of the tested optimizations, and the most ideal optimization has the lowest potential energy. Potential energy is measured as the MMFF energy -- a lower potential energy indicates a more stable molecule.

Calculate the MMFF energy of the optimized structure:

Previous optimizations with different lanraplenib positions were also run, which were formed by changing the placement of temporary hydrogen bonds between IL-8 and lanraplenib:

Grab previous optimizations from the cloud:

Graph the previous optimizations:

Calculate the previous optimizations’ MMFF energy:

The “valley” docked optimization forms the lowest MMFF energy out of all previous calculations. This optimization appears to interact with GLU-2, LEU-3, supporting a portion of the study’s pharmacophore analysis, but it interacts primarily with ARG-24, ILE-38, and GLN-6. While the study’s docking was similarly aligned as previous optimizations show above, this new docking has a much lower MMFF energy comparatively. This altered residue interaction is supported by the current understanding of interleukin-8’s binding site. In a 2020 study by Boyles et al., interleukin-8’s primary binding site, with certain “hot spot” residues being the primary mediators for interaction was identified (Boyles et al., 2020). This binding region’s residues includes ILE-38, LEU-3, and ARG-4, which all interacted with the created optimization of lanraplenib. The optimized interaction of lanraplenib shows an increased affinity to the binding site, rather than the previously identified residues. Lanraplenib’s interaction with the binding site of IL-8 points to a possible inhibitory capability.

Conclusion

Molecular and quantum mechanics techniques for multiple molecule energy optimization for both amino acid-small molecule interaction and protein-small molecule interaction was shown. From the Hartree-Fock and Density Functional Theory energy optimization computations of stachydrine and glutamic acid, 2 hydrogen bonds were shown to form. A method for successful small-molecule binding to proteins in Wolfram was created using MoleculeComplex[] and MoleculeManipulate3D[], with a streamlined expansion into residue and small molecule analysis. A Molecule object-based residue labelling system was also created to visualize residues and corresponding atomic indices. Modeling optimized interleukin-8 and lanraplenib binding using the Merck Molecular Force Field (MMFF) showed unique differences in binding site prediction compared to a pharmacophore model, indicates lanraplenib’s possibility as an interleukin-8 inhibitor.

Future Directions

Lanraplenib’s binding to interleukin-8’s identified binding site and particular “hot-spot” residues indicates more detailed analysis and a quantum mechanics approach to understand the chemical binding is needed. Different proteins and small molecules can be analyzed as well to find binding sites. Further study into protein formal charge effects when modeling in Wolfram is needed to better understand the importing and conversion of proteins. A way to change a ligand’s coordinates in relation to the protein to test multiple binding sites can be implemented to develop a systemic binding capability. From further computations, an electron density surface visualization can potentially be created to better visualize binding sites.

Protein-protein interaction can also be modeled, which can be used to confirm or discover the binding sites of interleukins, interleukin receptor complexes, and other proteins. Protein-protein interaction can also be modeled with a bound small-molecule to test inhibition and other binding properties.

Acknowledgements

I would like to thank my mentor, Dr. Jason Sonnenberg for their expert guidance and advising for this project, Dr. Jason Biggs for his work on the QuantumChemistry paclet and for allowing us to use it for this project, Dr. Robert Nachbar for his development of the MoleculeComplex paclet and technical advice, and the teaching assistants Ameya Srinivasan, Anne Shuai, Annie Zhu, Austin Zhang, Ian Gu, and William Choi-Kim for their continued help throughout the project.