Molecular docking software serves the purpose of assessing and predicting molecular interactions and recognition. It accomplishes this by revealing potential binding sites and modes at a structural level, as well as by calculating the binding affinity of these interactions. This application is commonly referred to as ligand-protein docking. Molecular docking is frequently employed to predict how small molecules interact with larger target molecules.

The receptor and ligand screening discussion can be found in the Project-Engineering section. After research and screening, the ideal ligand candidate was Honokiol. CB1 was deemed as the ideal receptor for the drug Honokiol, as AM-6538 is a strong CB1 antagonist, which would enable controlled release of these compounds from CB1.

We performed Molecular Docking Simulations to answer these questions:

1. Does Honokiol bind with CB1 to create a stable complex?
2. When AM-6538 is docked with CB1, does it bind at the same site as Honokiol?

We used AutoDock4 to perform these Molecular Docking simulations and PyMOL to visualize the complexes.

Receptor Preparation:

We used the 5tgz pdb, which is a crystal structure of the Human Cannabinoid Receptor CB1. The excess domain in the CB1 sequence and other small molecules in the crystal structure were removed using PyMOL. The domain was identified as the separate floating lobe opposite the free end of the primary group. The resultant is a coherent cylindrical structure functioning as the receptor, which is easier to engineer and fuse with other components. AutoDock was used to add polar hydrogens and repair the receptor residues with missing atoms.

Docking results between CB1 and Honokiol:

The highest binding strength was found to be -8.72 kcal/mol, followed by -8.54 kcal/mol. Honokiol was found to bind in a cylindrical opening in CB1, thus forming a secure binding. The region of binding for the conformations with the best score was common and can be highlighted from Figure 1.




Docking results between CB1 and AM-6538:

The complex achieved its best binding energy scores at -14.6 kcal/mol and -14.27 kcal/mol. The binding site for the same is depicted in Figure 2. Notably, we observed that AM-6538 occupies the identical binding site within CB1 as Honokiol. Consequently, based on its stronger binding energy and the common binding site, we can infer that AM-6538 functions as an antagonist, displacing Honokiol from the CB1 pocket.

To strengthen our inference, we conducted docking experiments with the CB1-Honokiol complex and the antagonist AM-6538. The outcomes revealed binding energy scores of -11.36 kcal/mol and -10.35 kcal/mol, indicating highly favorable competitive binding.

Molecular dynamics simulations allow us to understand the dynamic behavior of molecules and atoms. They play a pivotal role in fields like protein folding, and drug discovery, offering a window into the microscopic world of molecular interactions and motions.

We performed molecular dynamics simulations of the CB1-AM6538 complex in water to determine its stability. The system was simulated using the GROMACS MD engine and the AMBER99SB force field and explicit water(TIP3P water model). The simulation temperature was 313K and 0.15 Molar NaCl (with neutralizing ions) was used and Antechamber was used to create ligand topologies. Trajectory length was 100 nanoseconds.



To evaluate the simulated complex stability, Root mean square deviation (RMSD), a measure of protein structural drift, is often used. It indicates the sum of all atomic deviations of the conformation at a certain time and the target conformation. Each residue is compared to its position to a reference frame, and its RMSD is plotted for every frame. The RMSD is given by:

where, m_i is the mass of atom i, x_i is the coordinate vector of the target atom i, y_i is the coordinate vector of the reference atom i, and M is the total mass.

Over the 100 ns timeframe, the complex consistently maintains RMSD values below the 2Å threshold, a widely accepted criterion for molecular stability in molecular dynamics studies.

Mathematical modeling plays a pivotal role in advancing our understanding of tumors and their treatment. By leveraging mathematical principles and computational tools, these models uncover insights into tumor dynamics, growth patterns, and responses to treatment and also enable the optimization of personalized treatment strategies.

Mathematics can be used to guide laboratory research by anticipating the properties and behavior of components and systems. So, in order to forecast how certain components of our design would behave, we utilized a tumor growth mathematical framework. Our model takes inspiration from an in vivo experiment^[4] that includes a mixture of target antigen-positive and target antigen-negative tumor cells.

Our model involves the progression of four separate cellular groups throughout the treatment period. In particular, we analyze the entire tumor population as comprising two subgroups labeled as T_s and T_r, representing the antigen-positive and antigen-negative tumor cells present in the solid tumor, respectively. The other two cell populations include the CAR-T cells (C) and the Enhanced Immune Response cells (E).

The term enhanced immune response refers to the phenomenon where CAR-T cells trigger CD8 T cells, innate immune cells (such as macrophages, neutrophils, or natural killer (NK) cells), or even antigen-presenting cells (APCs) involved in antigen cross-presentation, allowing CAR T-cells to exert their therapeutic impact beyond the targeted antigen-expressing cells. This enhanced immune response enables CAR-T cells not only target cancer cells expressing a specific antigen but also stimulate cytotoxic activity against nearby cancer cells that lack the targeted antigen expression.

The equations:

1. \[ \frac{\mathrm{dT_{s}} }{\mathrm{d} t}=r_{1}T_{s}(1-\frac{T_{s}+T_{r}}{K_{1}})- D_{E}T_{s}-D_{C}T_{s}\]
2. \[ \frac{\mathrm{dT_{r}} }{\mathrm{d} t}=r_{2}T_{r}(1-\frac{T_{s}+T_{r}}{K_{1}})- D_{R}T_{r}\]
3. \[ \frac{\mathrm{dC} }{\mathrm{d} t}=I_{C}-\gamma_{C}C-\mu _{C}log(\frac{E+C}{K_{2}})\frac{D_{C}^{2}}{k+D_{C}^{2}}C - \omega_{C}CT{s}\]
4. \[ \frac{\mathrm{dE} }{\mathrm{d} t}=e-\gamma_{E}E-\mu _{E}log(\frac{E+C}{K_{2}})\frac{D_{E}^{2}}{k+D_{E}^{2}}E - \omega_{E}E(T{s}+T_{r})\]

where

\[D_{E}=d_{E}\frac{(E/T_{s})^l}{s+(E/T_{s})^l}\] \[D_{C}=d_{C}\frac{(C/T_{s})^l}{s+(C/T_{s})^l}\] \[D_{R}=d_{E}\frac{(E/T_{r})^l}{s+(E/T_{r})^l}\]

The first two equations (1) and (2) exhibit the logistic growth of tumor populations T_s and T_r, each with its own growth rate. The killing of the antigen-positive tumor population is modelled through two ratio-dependent cell-kill terms denoted by D_E and D_E. Those terms were originally proposed to capture the interaction between tumor cells and immunological components, including CD8+ T cells and natural killer cells. The effect on antigen-negative tumour population, T_r follows the same ratio-dependent cell-kill term denoted by D_R ^[5].

In the equations (3) and (4), the parameters γ_Cand γ_E represent the natural death rates of CAR T-cells and enhanced response cells, respectively. The parameter ω_C represents the rate of exhaustion of CAR T-cells resulting from exposure to the antigen-positive tumor population, while ω_E represents the rate of exhaustion of enhanced response cells due to their interaction with both tumor cell populations. We assume that the recruitment of CAR T-cells and enhanced response cells follows the Michaelis-Menten dynamics^[2]. The recruitment is represented by the terms D_C²/(k + D_C²) and D_E²/(k + D_E²), respectively. The recruitment of CAR-T cells and enhanced response cells are scaled by the immune cell competition term log((E+C)/K₂), which limits the proliferation of immune cell components beyond the carrying capacity K₂, with their maximum recruitment being denoted by μ_C and μ_E, respectively^[3].

CAR T-cell injection is modelled as a one-time increase in the amount of CAR T-cells on the days of injection. In equation (4), I_c represents the amount of CAR T-cell reaching the tumor site. In our model, we assume that 10⁶ cells can be approximated by 1mm³ of cells. Initial conditions are set as T_s(0) = 50p, T_r(0) = 50(1−p), C(0) = 0 and E(0) = 0.05 where p is the percentage of the antigen-positive tumor cell population. All units are given in mm³.

Parameters

The parameters used in the above set of equations along with their description and numerical values are tabulated below. In the absence of literature on the parameters, we have made educated guesses about some parameter values.

Paramter	Description	Value	Units
r1	Antigen postive cells proliferation rate	0.151	day-1
r2	Antigen negative cells proliferation rate	0.180	day-1
K1	Carrying capacities of tumor	5.9 x 103	mm3
l	exponent of tumor lysis by either C or E	1.291	unit-less
μC	max recruitment rate of CARTs by antigen-positive tumor lysis	0.6	day-1
dC	maximum killing rate of antigenpositive cells via CAR-T cells	0.27	day-1
s	steepness of fractional antigen-negative tumor kill by E cells	3.05 x 10-1	unit-less
γC	CAR-T's death rate	2.93 x 10-2	day-1
k	steepness of CAR-T/Enhanced cell recruitment	2.019 x 10-7	day-2
ωC	CAR-T exhaustion due to antigen-positive cells	3 x 10-5	mm-3 day-1
K2	immune cell carrying capacity	1.65 x 103	mm3
dE	maximum killing rate of antigen-positive/antigen-negative cells by enhanced response cells	0.27	day-1
μE	max recruitment rate of E by antigen-positive tumor lysis	0.82	day-1
e	base recruitment rate of E cells	5 x 10-2	mm3 day-1
γE	E cell death rate	2 x 10-2	day-1
ωE	E cell exhaustion due to antigen-positive/antigen-negative cells	3.42 x 10-6	mm3 day-1
u	Concentrion of drug	25	μM
a1	Fraction of antigen postive cells killed by drug	0.3	unit-less
a2	Fraction of antigen negative cells killed by drug	0.3	unit-less
a3	Fraction of C cells killed by drug	<a1,2	unit-less
a4	Fraction of E cells killed by drug	<a1,2	unit-less

Effect Of Drug:

We further improve our model by including the drug Honokiol and obtain its effect on the progression of the tumor populations of T_sand Tr. The amount of drug present at the tumour location at time t is designated as u(t). With the response curve in all cases being represented by an exponential F(u)= a(1- e^-ku) ^[5] and F(u) being the fraction of cells killed for a given dose of drug, u, at the tumour site. Since the details of the pharmacokinetics are unknown, we let k = 1 in these preliminary studies. We assume that the drug kills all types of cells but the kill rate varies for each type of cell.

The improved Equations are:

1. \[\frac{\mathrm{dT_{s}} }{\mathrm{d} t}=r_{1}T_{s}(1-\frac{T_{s}+T_{r}}{K_{1}})- D_{E}T_{s}-D_{C}T_{s} - a_{1}(1-e^{-u})T_{s}\]
2. \[\frac{\mathrm{dT_{r}} }{\mathrm{d} t}=r_{2}T_{r}(1-\frac{T_{s}+T_{r}}{K_{1}})- D_{R}T_{r}-a_{2}(1-e^{-u})T_{r}\]
3. \[\frac{\mathrm{dC} }{\mathrm{d} t}=I_{C}-\gamma_{C}C-\mu _{C}log(\frac{E+C}{K_{2}})\frac{D_{C}^{2}}{k+D_{C}^{2}}C - \omega_{C}CT{s}-a_{3}(1-e^{-u})C\]
4. \[\frac{\mathrm{dE} }{\mathrm{d} t}=e-\gamma_{E}E-\mu _{E}log(\frac{E+C}{K_{2}})\frac{D_{E}^{2}}{k+D_{E}^{2}}E - \omega_{E}E(T{s}+T_{r})-a_{4}(1-e^{-u})E\]

Our model centers on the tumor site subsequent to the allosteric displacement of Honokiol by AM-6538 from the 'seeker' T-cell, which occurs once these seeker T-cells have successfully reached the tumor site. CAR T cells can then be administered intravenously.

Analysis:

1. Tumor growth Model:
   We modeled the progression of tumor growth in a scenario where no treatment was applied. In this context, without the presence of CAR-T cells and the medication, the tumor mass exhibited continuous growth.
    
   Figure 5. Tumor progression without treatment.
2. CAR-T cell administration:
   This plot illustrates the impact of CAR-T cells and the resulting enhanced immune response on the tumor population. CAR-T cells demonstrate efficacy against antigen-positive tumor cells, but the antigen-negative cell population persists for a longer duration.
    
   Figure 6: Tumor progression with CAR-T cell administration.
3. Honokiol administration:
   Honokiol, when used in isolation, does not exhibit any significant positive impact on the tumor population. While a higher concentration of honokiol might lead to a reduction in tumor mass, it's important to note that the cytotoxic effects of honokiol impose an upper limit on the concentration that can be safely administered.
    
   Figure 7: Tumor progression with Honokiol administration.
4. CAR-T cells along with honokiol administration: The synergistic effect of CAR-T cells and honokiol results in a significant reduction in tumor mass within a shorter timeframe.
    
   Figure 8: Tumor progression with CAR-T cells along with honokiol.

Insights:

From our model, we derived the following insights: Independent CAR-T cells and honokiol treatment is not enough to control the tumor volume. The combined effect of CART and honokiol, however, leads to limiting the tumor mass and also significantly decreasing tumor volume in lesser time, showing more effectiveness. These findings can be analyzed to potentially reduce the required CAR-T cell concentration by regulating the honokiol concentration, thereby enhancing the cost-effectiveness of the treatment.

References:

1. Kara, E., Jackson, T.L., Jones, C., McGee II, R.L. and Sison, R., 2023. Mathematical Modeling Insights into Improving CAR T cell Therapy for Solid Tumors: Antigen Heterogeneity and Bystander Effects. arXiv preprint arXiv:2307.05606. Link
2. Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A. and Perelson, A.S., 1994. Nonlinear dynamics of immunogenic tumors: parameter estimation and global bifurcation analysis. Bulletin of mathematical biology, 56(2), pp.295-321.Link
3. Kimmel, G.J., Locke, F.L. and Altrock, P.M., 2019. Response to CAR T cell therapy can be explained by ecological cell dynamics and stochastic extinction events. bioRxiv, p.717074. Link
4. Klampatsa, A., Leibowitz, M.S., Sun, J., Liousia, M., Arguiri, E. and Albelda, S.M., 2020. Analysis and augmentation of the immunologic bystander effects of CAR T cell therapy in a syngeneic mouse cancer model. Molecular Therapy-Oncolytics, 18, pp.360-371. Link
5. De Pillis, L.G. and Radunskaya, A., 2001. A mathematical tumor model with immune resistance and drug therapy: an optimal control approach. Computational and Mathematical Methods in Medicine, 3(2), pp.79-100. Link
6. Qu, J., Mei, Q., Chen, L. and Zhou, J., 2021. Chimeric antigen receptor (CAR)-T-cell therapy in non-small-cell lung cancer (NSCLC): current status and future perspectives. Cancer immunology, immunotherapy, 70, pp.619-631. Link
7. Cheng, X., Wang, F., Qiao, Y., Chen, T., Fan, L., Shen, X., Yu, D., Huang, Y. and Wei, M., 2023. Honokiol inhibits interleukin-induced angiogenesis in the NSCLC microenvironment through the NF-κB signaling pathway. Chemico-Biological Interactions, 370, p.110295. Link