Overview

An important step in an experiment design is understanding how different reaction components interact and attaining an average time perspective of the experimental process.

Acknowledging this as a key point of our Model, we have addressed as many aspects of our project as possible.

We have integrated diverse elements including in-silico molecular design using Deep Learning algorithms, PDB files from X-Ray Crystallography, 3D structural analysis, Molecular Dynamics, Molecular Docking, and thorough Quantitative Timing Simulation Analysis.By doing so, we have successfully crafted a complex mathematical model, providing us with a comprehensive and holistic view of the behaviour of our biological system. This approach has enabled us to attain a profound understanding of the underlying mechanisms governing our system.

Figure 1 provides a detailed description of our system, with a focus on what is important to derive a mathematical model.

Figure 1: Schematic Representation of our System

Ordinary Differential Equations

Our mathematical model describes the molecular mechanisms of our experiment using Ordinary Differential Equations (ODEs).

Each chemical reaction step of our system is described by an ordinary differential equation.

The set of ODEs that mathematically characterize our system is the following:

The set of parameters that we chose from bibliography for our simulation is the following:

[1], [2], [3], [4], [5], [6], [7], [8],

System Analysis and Numerical Simulation

1^st step : NK cell transduction by lentivirus

Figure 2: NK cell transduction by lentivirus

The differential equation that describes the lentiviral transduction of NK cells, is the following:

We consider an initial NK cell concentration of 5 · 10⁶ cells per indicative volume unit (i.e., ml), which is 5 · 10⁶ cells/ml, and an initial lentiviral concentration of 1 MOI, which is 5.000.000 lentiviral particles for 5.000.000 NK cells (one lentivirus per NK cell).

For simplicity, we consider that an NK cell can be transduced by one lentivirus and that the NK cell lentiviral transduction success rate is 100% (the success rate depends on the experimental conditions).

If the success rate is, e.g., 20%, which means that 1.000.000 out of 5.000.000 NK cells got infected, then the variable values of the graph on the vertical axis are reduced to 20% of their initial value at any time moment of the horizontal axis.

The values of the following numerical simulation diagram are given in lentiviral RNA molecules per indicative volume unit (i.e., ml) on the vertical axis and in hours on the horizontal axis.

The graph of the 1st equation’s numerical simulation is the following:

Figure 3: Numerical simulation of the 1^st ODE

This graph represents the concentration of RNA molecules in the expanded NK cell culture after transduction, versus time.

[lentiviral RNA inside NK cell ] is the concentration of the lentiviral RNA molecules in the expanded NK cell culture after transduction occurs.

[lentiviral RNA outside NK cell ] is the concentration of the lentiviral RNA molecules inserted inside the expanded NK cell culture before transduction occurs.

[lentiviral RNA inside NK cell ] is affected by the process of reverse transcription that follows transduction, and that is represented by the negative term -Κ₂ · [lentiviral RNA inside NK cell ].

In the beginning (i.e., the first hours) of transduction, the concentration of lentiviral RNA in the expanded NK cell culture after transduction increases at a high rate, reaches a maximum value, and then decreases to practically zero value in about 75-80 hours, due to the reverse transcription of lentiviral RNA within each NK cell. Reverse transcription starts immediately after the moment when the first RNA molecules enter the cytoplasm of the NK cells in the expanded NK cell culture, i.e., at the moment t=0.

2^nd step : Reverse transcription of lentiviral RNA that entered NK cell cytoplasm and formulation of double-stranded lentiviral DNA

Figure 4: Lentiviral RNA reverse transcription after NK cell cytoplasm entry

The differential equation that describes reverse transcription of lentiviral RNA that entered NK cell cytoplasm and formulation of lentiviral DNA is the following:

[DNA before nucleus integration] is the lentiviral DNA concentration, and the lentiviral DNA molecules are located into the NK cell cytoplasm. Lentiviral DNA is the product of lentiviral RNA reverse transcription inside the NK cell cytoplasm.

The values of the following numerical simulation diagram are given in lentiviral DNA molecules per indicative volume unit (i.e., ml) on the vertical axis, and in hours on the horizontal axis.

The graph of the 2^nd equation’s numerical simulation is the following:

Figure 5: numerical simulation of the 2^nd ODE

The 2^nd differential equation substantially describes the process of reverse transcription of lentiviral RNA to double-stranded lentiviral DNA after transfection of the NK cells.

The negative term -Κ₃ · [DNA before nucleus integration] represents the simultaneous execution of the following reaction step, which is the lentiviral DNA integration into the NK cell nucleus. Due to this negative term, we observe the decrease of [DNA before nucleus integration] in the numerical simulation in Figure 5.

Lentiviral transduction begins at t=0 with initial [lentivirus concentration] = 1 MOI (i.e., [lentiviral RNA outside NK cell] = 5 · 10⁶ RNA molecules/ml), and reverse transcription starts immediately after the appearance of the first RNA molecules inside the cytoplasm of the NK cells i.e. at the same time t=0, with initial [lentiviral RNA inside NK cell] = 0.

Although reverse transcription is much faster than transduction (K_{reverse transcription} >> K_transduction), reverse transcription is completed shortly after transduction because it takes place after the transduction process.

Therefore, in the graph of Figure 5, the completion time of lentiviral RNA reverse transcription (horizontal axis) is a little longer compared to the graph of Figure 3, which represents the completion time of NK cells' lentiviral transduction.

This can be graphically observed in the following side-by-side enlarged graphs of Figures 3 and 5.

Figure 6 and 7: Enlarged graph of Figure 3 and Figure 5

Furthermore, the difference in the maximum value on the vertical axis of the graph of Figure 3 compared to the graph of Figure 5 is due to the fact that reverse transcription is much faster compared to transduction, while the subsequent process of reverse transcription, which is the lentiviral DNA integration into the NK cell nucleus, is much slower than reverse transcription (this can also be proved by comparing the corresponding K constants arithmetic values).

From the point of lentiviral DNA integration into the nucleus of the NK cell onwards, the following reactions (transcription of the lentiviral DNA, which is integrated into the nucleus of the NK cell, and translation of the lentiviral mRNA to CAR receptor protein) take place rapidly.

3^rd step : Lentiviral DNA integration into NK cell nucleus

Figure 8: Lentiviral DNA integration into NK cell nucleus

The differential equation that describes lentiviral DNA integration into NK cell nucleus is the following:

[DNA before nucleus integration] is the lentiviral DNA concentration inside the NK cell cytoplasm, and the lentiviral DNA molecules are located into the NK cell cytoplasm. Lentiviral DNA is the product of lentiviral RNA reverse transcription inside the NK cell cytoplasm.

The values of the following numerical simulation diagram are given in lentiviral DNA molecules per indicative volume unit (i.e., ml) on the vertical axis, and in hours on the horizontal axis.

Figure 9: Numerical simulation of the 3^rd ODE

We consider that the graph of Figure 9 corresponds to 100% lentiviral DNA nuclear integration success (i.e., integration of 5 · 10^>6 lentiviral DNA molecules into the nucleus of 5 · 10⁶ NK cells)

A crucial point that we must take into consideration is that lentiviral DNA integration into the NK cell nucleus takes place after two chemical reactions, which are the lentiviral transduction of the NK cell and the lentiviral RNA reverse transcription.

Therefore, at the first moments right after t=0, [DNA before nucleus integration] has a very small value because this is the time when the two preceding chemical reactions (lentiviral transduction of the NK cell and lentiviral RNA reverse transcription) have just started executing.

So, the speed value of the reaction

is pretty small at the first moments right after t=0, compared to the speed value afterwards, which means that the slope of the graph in Figure 9 is relatively small at the first moments right after t=0.

As a result, we expect that in the following graph, which represents the process of lentiviral mRNA translation, that this slope is going to be a little smaller (because the lentiviral mRNA translation comes after the lentiviral DNA integration into the NK cell nucleus and the lentiviral DNA transcription), and as a result, the speed of the lentiviral mRNA translation will have a little smaller value compared to the lentiviral DNA nuclear integration speed right after t=0.

Another thing that we observe is that the graph of Figure 3 (lentiviral transduction) has a steeper slope than the graph of Figure 5 (reverse transcription of lentiviral RNA), and the graph of Figure 5 (reverse transcription of lentiviral RNA) has a steeper slope than the graph of Figure 9 (Lentiviral DNA integration into NK cell nucleus) until they reach their maximum value, because the corresponding chemical reaction is waiting for the products of the previous reaction in order to proceed.

Of course, there is a chance that the lentiviral DNA integration into the NK cell nucleus is unsuccessful, which means that the lentiviral DNA has entered the NK cell nucleus and exits afterwards.

The following differential equation describes this phenomenon:

4^th step : Lentiviral DNA escape from the NK cell nucleus

Figure 10: Lentiviral DNA escape from the NK cell nucleus

The differential equation that describes lentiviral DNA nuclear escape is the following:

The following numerical simulation diagram values are given in lentiviral DNA molecules per indicative volume unit (i.e., ml) on the vertical axis and in hours on the horizontal axis.

Figure 11: Numerical simulation of the 4^th ODE

The graph of Figure 11 describes the case where ten lentiviral DNA molecules entered the nucleus of ten NK cells but finally did not achieve nuclear integration and escaped afterwards. They escape after almost 8-10 hours.

The number of lentiviral DNA molecules that do not achieve nuclear integration and escape the nucleus afterwards varies according to the experimental conditions and the reaction rate K, and also the escape time of each lentiviral DNA molecule may differ.

5^th step : Transcription of the lentiviral DNA, which successfully integrated into the NK cell nucleus, production of lentiviral mRNA, and degradation of lentiviral mRNA

Figure 12: Transcription of the lentiviral DNA, which successfully integrated into the NK cell nucleus, production of lentiviral mRNA, and degradation of lentiviral mRNA

The differential equation that describes transcription of the lentiviral DNA (which successfully integrated into the NK cell nucleus), the production of CAR gene mRNA molecules, and degradation of CAR gene mRNA molecules, is the following:

The values of the following numerical simulation diagram are given in lentiviral mRNA molecules per indicative volume unit (i.e., ml) on the vertical axis, and in hours on the horizontal axis.

Figure 13: Numerical simulation of the 5^th ODE

The graph of Figure 13 describes the variation in the concentration of the CAR gene mRNA ([mRNA of CAR gene]) during the transcription process of the CAR gene, which integrated successfully into the NK cell nucleus, and during the process of CAR gene mRNA degradation after transcription.

The negative term -Κ_{mRNA_deg} · [mRNA of CAR gene] represents the simultaneous execution of the following reaction step, which is mRNA degradation (mRNA survives for a finite period of time in the cytoplasm). Due to this negative term, we observe the decrease of [mRNA of CAR gene] in the numerical simulation in Figure 13.

Even though transcription is much faster than the process of lentiviral DNA integration into the nucleus of the NK cell (K_transcription >> K_nuclear_entry), the completion time of transcription is slightly longer compared to the completion time of lentiviral DNA integration into the nucleus of the NK cell (practically, they are equal) because transcription follows the process of lentiviral DNA integration into the nucleus of the NK cell.

6^th step : CAR gene mRNA translation to CAR receptor protein, and CAR receptor protein degradation

Figure 14: CAR gene mRNA translation to CAR receptor protein, and CAR receptor protein degradation

The differential equation that describes CAR gene mRNA translation to CAR receptor protein and CAR receptor protein degradation is the following:

The values of the following numerical simulation diagram are given in CAR receptor molecules per indicative volume unit (i.e., ml) on the vertical axis, and in hours on the horizontal axis.

Figure 15: Numerical simulation of the 6^th ODE

The CAR receptor degradation rate is very low (as shown by the small value of K_{protein_deg}), and therefore, practically during the time interval needed for CAR receptor to neutralize the cancer cells, its concentration, which is [CAR receptor] , is slightly reduced due to degradation.

As mentioned before (at the stage of DNA integration into the nucleus of the NK cell), we expect the graph slope close to t=0 (at the very first moments right after t=0) to be a little bit smaller. The reason is that lentiviral mRNA translation follows lentiviral DNA transcription, and the speed of lentiviral mRNA translation will be a little lower than the lentiviral DNA integration speed at the first moments right after t=0.

Therefore, we observe that CAR receptor protein production is completed in about 75-80 hours, which is about 3 days.

The choice of K constants derived from experimental confirmation can approximate the reaction completion time more accurately.

In the following steps, we model the CAR receptor binding reaction with the mesothelin antigen found on the surface of pancreatic cancer cells.

7^th step : CAR receptor binding with mesothelin antigen

Figure 16: CAR receptor binding with mesothelin antigen

The differential equation that describes CAR receptor binding with mesothelin antigen is the following:

We consider that the initial CAR receptor protein concentration is equal to 10^-6 molarity ([CAR receptor]=10^-6 M), and initial mesothelin antigen concentration is equal to 10^-6 molarity ([mesothelin antigen]=10^-6 M)

The values of the following numerical simulation graphs are given in CAR receptor molarity (M), mesothelin antigen molarity (M), and CAR receptor-antigen complex molarity (M) on the vertical axis, and in seconds on the horizontal axis.

We observe that CAR receptor binding to the antigen is completed in about 10 seconds.

What should be mentioned here is that the higher the initial concentrations of the reacting CAR protein and antigen molecules ([CAR receptor], [mesothelin antigen]) are, the higher the reaction speed is and the shorter the reaction completion time is. This is because the number of collisions of reacting molecules per unit time increases, and therefore, the number of effective collisions per unit time (i.e., colliding molecules with proper orientation) increases.

Furthermore, the reaction speed also depends on the value of the constant K₈.

So if the reaction has a high K₈ (even with low antigen concentration), it will make sure that our CAR-NK cells will bind to target antigen fast, even when the antigen concentration is low [6].

In addition, an essential point regarding the binding reaction between CAR receptor and antigen is that for concentration values of [mesothelin antigen] = 10^-7 M and below, the binding reaction time is stabilized to 1 second.

This can be observed from the following graphs during the numerical simulation of the 7^th equation:

So, the optimal rate of mesothelin antigen concentration and CAR receptor concentration is

8^th step : Cancer cell neutralization by CAR-NK cell

The differential equation that describes cancer cell neutralization by CAR-NK cell is the following:

We consider that the initial concentration of CAR-NK - cancer cell complex is equal to 10^-6 molarity (M)

The values of the following numerical simulation graph are given in CAR NK-cancer cell complex molarity (M) on the vertical axis, and in hours on the horizontal axis.

So, we observe that the death of pancreatic cancer cells that have mesothelin on their surface by CAR-NK cells takes place in about 5 hours.

With a suitable choice of a K constant for this reaction derived from experimental data, the death rate of pancreatic cancer cells can be approximated accurately.

In Silico Chimeric Antigen Receptor Design

Our initial design was conceived based on the research conducted by Ye Li et al. [9], who systematically assessed ten distinct chimeric antigen receptors (CARs) for their stability and their capacity to selectively and efficiently target cancer cells. Among their evaluated models, the most promising CAR construct featured the following components:

1.scFv : A single-chain variable fragment specific to the mesothelin antigen.

2.NKG2D : An intramembrane domain.

3.2B4-CD3ζ : An intracellular domain responsible for signal transduction upon mesothelin antigen recognition.

Armed with this knowledge, our research aimed to computationally generate this CAR structure "in silico". Initially, we attempted to locate the complete structure in the RCSB Protein Data Bank but were unsuccessful. Consequently, we created the full tertiary structure de novo, leveraging the computational tools available within the Robetta platform [10].

Robetta offers rapid and precise deep learning-based methodologies, such as RoseTTAFold and TrRosetta, along with an interactive submission interface that facilitates custom sequence alignments for homology modeling, constraints, local fragments, and other parameters. It supports the modeling of multi-chain complexes using RoseTTAFold (with paired MSA provided by the user) or comparative modeling (CM) and provides options for large-scale sampling.

We retrieved the amino acid sequences of each individual CAR component from the NCBI database. We assembled them into a single input for the Robetta web server, resulting in the generation of the complete tertiary structure of the CAR. To ensure the accuracy of our results, we cross-validated each constituent component with its corresponding entry in the Protein Data Bank.

Furthermore, to closely replicate the in vivo behavior of our CAR, we aimed to create a model of the cell membrane surrounding it. We derived inspiration from the work of Irina D. Pogozheva et al. [11] , particularly their eukaryotic plasma membrane model. The construction of the membrane was executed using the CHARMM-GUI [12] [13] [14] [15] [16] [17] [18] [19], employing the following lipid composition (where '#' represents numerical values):

In this manner, we successfully created the Chimeric Antigen Receptor specific for Mesothelin.

Expanding upon our model, we investigated the interaction between the CAR and its target, the Mesothelin antigen. We obtained the full-length Mesothelin structure from the Protein Data Bank using the PDB code '7UED' [20], extracting only the antigen and removing the antibody structure.

Molecular Dynamics Simulations

Subsequently, we employed GROMACS [21] [22], a versatile Molecular Dynamics software for our Molecular Dynamics simulations. The computational procedures employed in our study adhered to rigorous scientific standards and followed a well-established protocol. We present a comprehensive description of the critical simulation parameters and methodologies below:

Force Field Selection: We opted for the CHARMM36m force field, as detailed in a seminal publication by Huang and MacKerell [23]. To represent the protein, this force field was meticulously chosen for its capacity to accurately capture the intricate interactions governing protein behavior. Furthermore, the TIP3P water model, originally developed by Cornell et al. [24], was employed to characterize the water molecules within the simulation environment.

Simulation Box: We constructed a cubic simulation box with the molecule of interest positioned at its geometric center. A distance of 2.0 nanometers (nm) was maintained between the molecule and the box boundaries to minimize boundary effects and provide an adequate buffer. To neutralize the system's overall charge, sodium ions (Na+) and chloride ions (Cl-) were judiciously introduced into the simulation.

Energy Minimization: The optimization of energy within the system was achieved through an energy minimization process. This endeavor employed a steepest descent algorithm with a maximum force tolerance (emtol) set at 1000.0 kJ · mol^-1 · nm^-1. This step aimed to alleviate unfavorable interactions and attain a stable starting configuration for subsequent simulations.

Temperature Equilibration (NVT): The simulation was executed under canonical ensemble conditions (NVT), where the number of particles (N), volume (V), and temperature (T) were held constant. The temperature was meticulously regulated at 310.15 Kelvin (37.0°C) using the Berendsen thermostat. Coupling groups were designated for the relevant entities, namely CAR, MSLN (Mesothelin), and SOLV (solvent). This temperature control was achieved through a thermostat with a temperature time constant of 0.2 picoseconds (ps). To ensure a gradual and controlled equilibration process, three successive rounds of NVT equilibration were conducted, progressively reducing restraints on the molecules.

Pressure Equilibration (NPT): The isothermal-isobaric ensemble (NPT) was adopted to maintain a constant pressure of 1 bar. This was accomplished using a Parrinello-Rahman barostat, ensuring that the simulation system remained at physiologically relevant pressure conditions.

In adhering to these meticulously selected parameters and methodologies, we were able to perform our simulations with scientific rigor, thereby ensuring the reliability and reproducibility of our findings. To facilitate the discrimination of distinct elements within our simulation, we established a comprehensive index: residues 0-287 were designated to represent the Mesothelin antigen, residues 288-3015 corresponded to the Chimeric Antigen Receptor (CAR) intricately associated with the cell membrane, and residues 3016-257310 were allocated to signify the surrounding solvent molecules. The results of our exhaustive 8-nanosecond simulation are meticulously presented here.

Evidently, throughout the entire simulation duration, the Mesothelin antigen consistently maintained its proximity to the CAR, thereby affirming the anticipated interaction that was initially hypothesized. Furthermore, a meticulous examination of the formation of Hydrogen Bonds (HB) during the simulation period between the CAR and the Mesothelin antigen revealed compelling insights. As the simulation progressed over time, a discernible augmentation in the number of formed Hydrogen Bonds was observed, transitioning from an initial count of 6600 to a subsequent tally of 7000, as visually depicted in the graph presented below. These newly established bonds between the receptor and the antigen distinctly substantiate the emergence of a meaningful interaction between the receptor and the C-terminal of the antigen.

Antibody Optimization

Binding between antibody and antigen is developed between antibody’s third complementarity-determining region (CDR3) on the heavy chain, and the antigen surface.

The most important step for developing effective therapeutic antibodies is to design CDRs that bind to the specific antigen [25]

In our antibody optimization study, we analyzed the complex of MORAb 15B6 antibody - Mesothelin antigen [26], which is represented in the PBD file ‘7u8c’ in RCSB [27] and in the SabDab Database [28].

The antibody optimization was conducted with the use of Diffab [25], a Diffusion-Based Generative Model for Antigen-Specific Antibody Design and Optimization.

In order to optimize an antibody, the CDR-H3 sequence and structure was first perturbed for t steps using the forward diffusion process.

Then, the sequences starting from the (T − t)-th step (t steps remaining) of the generative diffusion process were denoised , and then a set of 100 optimized antibodies was obtained.

The optimization process was conducted for t=64 steps, and the best optimized antibodies with the best (lowest) binding energy were those in the t=4 and t=8 optimization step [25].

The optimization steps t=16, t=32 and t=64 produced antibodies which had many changes in the CDR-H3 amino-acid sequence, which resulted in high energy, instability of the antibody-antigen complex and significant changes in the coordinates of the atoms of the CDR-H3 amino-acids, whereas, the t=2 step produced antibodies with no change in their CDR-H3 amino-acid sequence, only just some changes in the coordinates of the atoms of the CDR-H3 amino-acids which resulted in no significant binding energy changes.

In t=4 and t=8 optimization step, the CDR-H3 amino-acid sequence was very similar to the original one, and only some amino-acids in some of the produced antibodies finally changed, which was a positive factor about the binding energy of new complex

The following graphs show the number of amino-acids that changed in every optimized antibody in t=4 and t=8 optimization step.

The initial CDR-H3 amino-acid sequence of our antibody is ELGGY.

In t=4 and t=8 optimization step there were some CDR-H3 amino-acid sequence changes in some of the produced antibodies.

The following tables show these changes.

Afterwards, we conducted a binding energy study for the antibody-antigen complexes in t=4 and t=8 step.

Binding energy of an antibody-antigen complex is equal to the total energy of the complex minus the sum of the protein and ligand energy before binding.

In other words, it is described by the following equation:

and should have a negative value, because binding energy is the energy released into the environment by the creation of CAR protein-antigen binding bonds, i.e. it is the reduction of the total energy (thermal content) of the system, and is expressed as, for example, ΔH = -100 kcal/mol. ΔH is the enthalpy value change of the complex.

If a protein-ligand complex has a high value energy score, then it is not favorable, because it shows that it is unstable. This means that the higher the energy score is, the lower the stability is.

Therefore, we analyzed our produced complexes, and searched for those which had the lowest binding energy.

The software that the utilized for our energy analysis is PyRosetta [29], and the score function that we used had the following parameter weight values [30] :

First, we cleaned our produced complexes, so as to remove any water molecules from their PDB structure, and then we calculated their binding energy using our score function.

Our initial PDB file’s binding energy score was -62.182 kcal/mol.

The results that we obtained from our binding energy study for t=4 and t=8 step respectively are the following:

So, the best (lowest) binding energy using this score function was on the 25th PDB file produced by Diffab, which is -51,976 kcal/mol.

We also calculated the binding energy of the 25th PDB file using another score function, which is the latest Rosetta Energy Function, REF15 [30], , which has the following parameter weight values:

The binding energy score that we obtained from the REF15 energy function was -53,204 kcal/mol, which is better score compared to the previous score function.

This means that the energy function weight parameters play a significant role in our energy calculations, and different score functions can result in different energy scores.

However, in our case study, we conclude that our initial CAR (speicifically, the CDR-H3 part) has very good binding characteristics.

As a result, the computational energy calculations are a significant tool , but further wet lab experiments are urgent to verify our optimized antibodies. It remains unclear whether the generated antibodies can be produced in the wet lab and actually bind to the target. More efforts are needed to design a biologically effective antibody [25].

Molecular Docking Analysis

In addition to the comparative assessment of free energy between the initial CAR and the enhanced CAR design, we conducted a rigorous docking analysis involving the interaction with the Mesothelin antigen. This investigation was facilitated through the utilization of HDOCK [31], a web-based platform distinguished for its application of a hybrid strategy in protein-protein docking simulations. The outcomes yielded congruent findings with the energy-based comparison, as articulated below:

For the initial CAR design :

- Docking Score: -251.46

- Confidence Score: 0.8838

- Ligand RMSD: 103.07

For the optimized CAR design :

- Docking Score: -228.61

- Confidence Score: 0.8281

- Ligand RMSD: 116.65

It is imperative to underscore that higher docking scores, in conjunction with elevated confidence scores, signify a more robust binding affinity between the receptor and the ligand molecule. In light of these discerning results, it is unequivocally evident that the original CAR design exhibits superior binding characteristics, rendering it the most auspicious candidate poised to yield the most promising outcomes within this context.

It is of significance to acknowledge that the robust validation of these findings necessitates the execution of a comprehensive dynamic analysis concerning the binding dynamics between the enhanced CAR and the Mesothelin antigen, akin to the meticulous scrutiny applied to the initial CAR design. Regrettably, time restrictions precluded us from undertaking such elaborate experiments at this juncture. Furthermore, it is judicious scientific practice to extend the existing dynamics simulations, thereby fortifying the veracity and comprehensiveness of our research outcomes.

Future Steps

Our suggestion is the development of an extended mathematical model that will describe our experiment in-depth, in terms of atom and molecular scale (i.e., modeling the entrance of the lentivirus into the cell membrane of the NK cell in atom-level scale, and therefore extend our model to an in-depth mathematical representation of every molecular phenomenon that takes place in each experimental step). Also, the use of reaction rate constants K derived from wet lab experimental data would definetely enhance the precision of our computations and simulations. In other words, this would be a new way of mathematical modeling, including in-depth mathematical analysis of every molecular interaction step of our experiment, supposing that each molecular interaction is an independent mathematical sub-model itself. So, the experiment's total mathematical model would consist of several molecular mathematical sub-models. This extended model would have remarkable precision, would be very useful to check what could go wrong in each molecular step, and be an in-silico molecular testing and verification model for the wet lab experiment.

References

[1] Nguyen, T., Sha, S., Hong, M. S., Maloney, A. J., Barone, P. W., Neufeld, C., Wolfrum, J., Springs, S. L., Sinskey, A. J., & Braatz, R. D. (2021). Mechanistic model for production of recombinant adeno-associated virus via triple transfection of HEK293 cells. Molecular Therapy - Methods & Clinical Development, 21, 642–655.

2] https://assets.thermofisher.com/TFS-Assets/LSG/brochures/superscript-alternative-amv-white-paper.pdf

[3] Team:Northwestern/Model - 2018.igem.org. (n.d.). https://2018.igem.org/Team:Northwestern/Model

[4] Median mRNA decay constant - Budding yeast Saccharomyces ce - BNID 105510. (n.d.). https://bionumbers.hms.harvard.edu/bionumber.aspx?s=n&v=5&id=105510

[5] Search BioNumbers - The Database of Useful Biological Numbers. (n.d.). https://bionumbers.hms.harvard.edu/

[6] Mao, R., Kong, W., & He, Y. (2022). The affinity of antigen-binding domain on the antitumor efficacy of CAR T cells: Moderate is better. Frontiers in Immunology, 13. 

[7] Liu, J., Ma, C., Zhang, Z., & Chen, W. (2021). A computational model of CAR T-cell immunotherapy predicts leukemia patient responses at remission, resistance, and relapse. medRxiv (Cold Spring Harbor Laboratory).

[8] Sahoo, P., Yang, X., Abler, D., Maestrini, D., Adhikarla, V., Frankhouser, D., Cho, H., Machuca, V., Wang, D., Barish, M. E., Gutova, M., Branciamore, S., Brown, C. E., & Rockne, R. C. (2020). Mathematical deconvolution of CAR T-cell proliferation and exhaustion from real-time killing assay data. Journal of the Royal Society Interface, 17(162), 20190734.

[9] Li, Y., Hermanson, D., Moriarity, B. S., & Kaufman, D. S. (2018). Human iPSC-Derived Natural Killer Cells Engineered with Chimeric Antigen Receptors Enhance Anti-tumor Activity. Cell Stem Cell, 23(2), 181-192.e5.

[10] Robetta - Baker Lab

[11] Pogozheva, I. D., Armstrong, G. A., Kong, L., Hartnagel, T. J., Carpino, C. A., Gee, S. E., Picarello, D. M., Rubin, A., Lee, J., Park, S., Lomize, A. L., & Im, W. (2022). Comparative molecular dynamics simulation studies of realistic eukaryotic, prokaryotic, and archaeal membranes. Journal of Chemical Information and Modeling, 62(4), 1036–1051.

[12] Jo, S., Kim, T., Iyer, V., & Im, W. (2008). CHARMM-GUI: A web-based graphical user interface for CHARMM. Journal of Computational Chemistry, 29(11), 1859–1865.

[13] Brooks, B. R., Brooks, C. L., MacKerell, A. D., Nilsson, L., Petrella, R. J., Roux, B., Won, Y., Archontis, G., Bartels, C., Boresch, S., Caflisch, A., Caves, L. S. D., Cui, Q., Dinner, A. R., Feig, M., Fischer, S., Gao, J., Hodošček, M., Im, W., . . . Karplus, M. (2009). CHARMM: The biomolecular simulation program. Journal of Computational Chemistry, 30(10), 1545–1614.

[14] Lee, J., Cheng, X., Swails, J., Yeom, M. S., Eastman, P., Lemkul, J. A., Wei, S., Buckner, J., Jeong, J. C., Qi, Y., Jo, S., Pande, V. S., Case, D. A., Brooks, C. L., MacKerell, A. D., Klauda, J. B., & Im, W. (2015). CHARMM-GUI Input Generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 Additive Force Field. Journal of Chemical Theory and Computation, 12(1), 405–413.

[15] Wu, E. L., Cheng, X., Jo, S., Rui, H., Song, K., Dávila-Contreras, E. M., Qi, Y., Lee, J., Monje-Galvan, V., Venable, R. M., Klauda, J. B., & Im, W. (2014). CHARMM-GUIMembrane Buildertoward realistic biological membrane simulations. Journal of Computational Chemistry, 35(27), 1997–2004.

[16] Jo, S., Lim, J. B., Klauda, J. B., & Im, W. (2009). CHARMM-GUI membrane builder for mixed bilayers and its application to yeast membranes. Biophysical Journal, 97(1), 50–58.

[17] Jo, S., Kim, T., & Im, W. (2007). Automated builder and database of Protein/Membrane Complexes for molecular Dynamics simulations. PLOS ONE, 2(9), e880.

[18] Lee, J., Patel, D. S., Ståhle, J., Park, S., Kern, N. R., Kim, S., Lee, J., Cheng, X., Valvano, M. A., Holst, O., Knirel, Y. A., Qi, Y., Jo, S., Klauda, J. B., Widmalm, G., & Im, W. (2018). CHARMM-GUI Membrane Builder for Complex Biological Membrane Simulations with Glycolipids and Lipoglycans. Journal of Chemical Theory and Computation, 15(1), 775–786.

[19] Park, S., Choi, Y. K., Kim, S., Lee, J., & Im, W. (2021). CHARMM-GUI Membrane Builder for Lipid Nanoparticles with Ionizable Cationic Lipids and PEGylated Lipids. Journal of Chemical Information and Modeling, 61(10), 5192–5202.

[20] https://www.rcsb.org/structure/7UED

[21] Van Der Spoel, D., Lindahl, E., Hess, B., Groenhof, G., Mark, A. E., & Berendsen, H. J. C. (2005). GROMACS: Fast, flexible, and free. Journal of Computational Chemistry, 26(16), 1701–1718.

[22] Abraham, M., Murtola, T. J., Schulz, R., Páll, S., Smith, J. C., Hess, B., & Lindahl, E. (2015). GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. SoftwareX, 1–2, 19–25.

[23] Huang, J., Rauscher, S., Nawrocki, G., Ran, T., Feig, M., De Groot, B. L., Grubmüller, H., & MacKerell, A. D. (2016). CHARMM36m: an improved force field for folded and intrinsically disordered proteins. Nature Methods, 14(1), 71–73.

[24] Cornell, W. D., Cieplak, P., Bayly, C. I., Gould, I. R., Merz, K. M., Ferguson, D. M., Spellmeyer, D. C., Fox, T., Caldwell, J. W., & Kollman, P. A. (1995). A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society, 117(19), 5179–5197.

[25] Luo, S., Su, Y., Peng, X., Wang, S., Peng, J., & Ma, J. (2022). Antigen-Specific Antibody Design and Optimization with Diffusion-Based Generative Models for Protein Structures. bioRxiv (Cold Spring Harbor Laboratory).

[26] Liu, X., Onda, M., Watson, N., Hassan, R., Ho, M., Bera, T. K., Jiang, W., Chakraborty, A., Beers, R., Zhou, Q., Shajahan, A., Azadi, P., Zhan, J., Xia, D., & Pastan, I. (2022b). Highly active CAR T cells that bind to a juxtamembrane region of mesothelin and are not blocked by shed mesothelin. Proceedings of the National Academy of Sciences of the United States of America, 119(19).

[27] https://doi.org/10.2210/pdb7U8C/pdb

[28] SABDAB: The Structural Antibody Database.

[29] Chaudhury, S., Lyskov, S., & Gray, J. J. (2010). PyRosetta: a script-based interface for implementing molecular modeling algorithms using Rosetta. Bioinformatics, 26(5), 689–691.

[30] Alford, R. F., Leaver‐Fay, A., Jeliazkov, J. R., O’Meara, M. J., DiMaio, F., Park, H., Shapovalov, M. V., Renfrew, P. D., Mulligan, V. K., Kappel, K., Labonte, J. W., Pacella, M. S., Bonneau, R., Bradley, P., Dunbrack, R. L., Das, R., Baker, D., Kuhlman, B., Kortemme, T., & Gray, J. J. (2017). The Rosetta All-Atom Energy Function for macromolecular modeling and design. Journal of Chemical Theory and Computation, 13(6), 3031–3048.

[31] Yan, Y., Tao, H., He, J., & Huang, S. (2020). The HDOCK server for integrated protein–protein docking. Nature Protocols, 15(5), 1829–1852.