Our Targets

Transferrin receptor

Melanotransferrin receptor

CAIX

α-Hemolysin

HBV surface antigen

Transferrin receptor binds to transferrin. In our simulation, we used the ectodomain of the transferrin receptor (in wheat color) with PDB ID: 1CX8¹.

Melanotransferrin receptor (LRP), which belongs to the low-density lipoprotein receptor protein family, binds to melanotransferrin. In our simulation, we used a binding domain of LRP (in firebrick color) with PDB ID: 2FYJ^2,3.

We used the HBV surface antigen (in chocolate color) with PDB ID: 7TUK⁴.

We also included alpha hemolysin (in pink color) with PDB ID: 7AHL⁵.

CAIX (in orange color) with PDB ID: 6FE2 was also used in our simulation⁶.

Structure and Docking prediction

Introduction

Structure prediction and docking simulation visualize the shape of the antibody and let us understand the antibody better,which is beneficial to further develop antibody modifications. Since the shape of a protein is informative for its function, the prediction of the designed antibody structure could suggest the optimal antibody construct candidates.

The goal of the docking simulation is to evaluate designed constructs in silico, use bioinformatics tools to simulate the actual performance of the antibody-target binding in order to select the highly functional construct to compare with the experimental data.

Approach: Step1: Design Antibodies with different constructs with the ability of binding to 1 or 2 targets

Step2: Predict the antibody structure with Alphafold2 using Colabfold⁷ and the variable region (Heavy/Light chain) of the antibody using SAbPred⁸.

Step3: Perform docking simulation of SAbPred predicted structure with targets using Hdock⁹ to simulate the docking for every target-variable region pairs.

Step4: Use PyMol (PyMOL Molecular Graphics System, Version 2.0 Schrödinger, LLC.) to superimpose the variable region on the predicted antibody structure to illustrate the structure-target docking

List of Constructs

Docking_score (Hdock): The docking scores are calculated by knowledge-based iterative scoring function ITScorePP or ITScorePR from Huanglab. A more negative docking score means a more possible binding model, but the score should not be treated as the true binding affinity of two molecules because it has not been calibrated to the experimental data. Protein-protein complex in PDB normally have a docking score of around -200 or better

Ligand RMSD (Hdock): The ligand RMSDs are calculated by comparing the ligands in the docking models with the input or modeled structures.

Alignment RMSD (PyMOL): Alignment RMSD is computed by a two-step approach for aligning structures: first, it performs a sequence alignment, and then it minimizes the Root Mean Square Deviation (RMSD) between the aligned residues.

C

Target: alpha hemolysin.

Docking_score (alpha hemolysin)	=	-250.37
Ligand RMSD (alpha hemolysin)	=	31.37
Alignment RMSD (alpha hemolysin)	=	0.795

D

Target: transferrin receptor and alpha hemolysin. Transferrin single-chain variable fragment attached N-terminally to anti-toxin A (hemolsyin) light and heavy chains of C.

Docking_score (alpha hemolysin)	=	-265.74
Ligand RMSD (alpha hemolysin)	=	56.5
Alignment RMSD (alpha hemolysin)	=	0.476
docking_score (transferrin)	=	-273.95
Ligand RMSD (transferrin)	=	31.37
Alignment RMSD (transferrin)	=	0.795

E

Target: transferrin receptor and alpha hemolysin. Transferrin single-chain variable fragment attached C-terminally to anti-toxin A (hemolsyin) light and heavy chains of C.

Docking_score (alpha hemolysin)	=	-248.49
Ligand RMSD (alpha hemolysin)	=	39.63
Alignment RMSD (alpha hemolysin)	=	0.541
docking_score (transferrin)	=	-283.84
Ligand RMSD (transferrin)	=	53.07
Alignment RMSD (transferrin)	=	0.922

F

Target: alpha hemolysin. Like construct C but light and heavy chain are connected by a linker.

Docking_score (alpha hemolysin)	=	-248.49
Ligand RMSD (alpha hemolysin)	=	39.63
Alignment RMSD (alpha hemolysin)	=	0.493

G

Target: transferrin receptor and alpha hemolysin. Like construct D but light and heavy chain are connected by a linker.

Docking_score (alpha hemolysin)	=	-250.37
Ligand RMSD (alpha hemolysin)	=	31.37
Alignment RMSD (alpha hemolysin)	=	0.555
docking_score (transferrin)	=	-287.47
Ligand RMSD (transferrin)	=	51.11
Alignment RMSD (transferrin)	=	0.890

H

Target: transferrin receptor and alpha hemolysin. Like construct E but light and heavy chain are connected by a linker.

Docking_score (alpha hemolysin)	=	-250.37
Ligand RMSD (alpha hemolysin)	=	31.37
Alignment RMSD (alpha hemolysin)	=	0.511
docking_score (transferrin)	=	-277.93
Ligand RMSD (transferrin)	=	48.56
Alignment RMSD (transferrin)	=	0.888

I

Target: alpha hemolysin. Like construct F but constant light chain is removed.

Docking_score (alpha hemolysin)	=	-250.37
Ligand RMSD (alpha hemolysin)	=	31.37
Alignment RMSD (alpha hemolysin)	=	000

J

Target: alpha hemolysin. Anti-hTfR single-chain variable fragment attached to anti-toxin A (hemolsyin) single-chain variable fragment fused to the heavy chain constant region

Docking_score (alpha hemolysin)	=	-241.84
Ligand RMSD (alpha hemolysin)	=	38.10
Alignment RMSD (alpha hemolysin)	=	0.00
Docking_score (hTfR)	=	-287.47
Ligand RMSD (hTfR)	=	51.11
Alignment RMSD (hTfR)	=	0.00

K

Target: melanotransferrin receptor and alpha hemolysin. Melanotransferrin peptide attached N-terminally to anti-toxin A (hemolsyin) light and heavy chains of C.

Docking_score (alpha hemolysin)	=	-241.84
Ligand RMSD (alpha hemolysin)	=	38.10
Alignment RMSD (alpha hemolysin)	=	0.00
Docking_score (melanotransferrin)	=	-143.38
Ligand RMSD (melanotransferrin)	=	88.34
Alignment RMSD (melanotransferrin)	=	0.00

L

Target: melanotransferrin receptor, alpha hemolysin. Melanotransferrin peptide attached C-terminally to anti-toxin A (hemolsyin) light and heavy chains of C.

Docking_score (alpha hemolysin)	=	-248.49
Ligand RMSD (alpha hemolysin)	=	39.63
Alignment RMSD (alpha hemolysin)	=	0.00
Docking_score (melanotransferrin)	=	-143.38
Ligand RMSD (melanotransferrin)	=	88.34
Alignment RMSD (melanotransferrin)	=	0.00

M

Target: melanotransferrin receptor, alpha hemolysin. Like K, but light chain and heavy chain are connected by a linker.

Docking_score (alpha hemolysin)	=	-248.49
Ligand RMSD (alpha hemolysin)	=	39.63
Alignment RMSD (alpha hemolysin)	=	0.00
Docking_score (melanotransferrin)	=	-143.38
Ligand RMSD (melanotransferrin)	=	88.34
Alignment RMSD (melanotransferrin)	=	0.00

R

Target: CAIX.

Docking_score (CAIX)	=	-291.98
Ligand RMSD (CAIX)	=	71.83
Alignment RMSD (CAIX)	=	0.00

W

Target: HBsAg.

Docking_score (HBsAg)	=	-355.62
Ligand RMSD (HBsAg)	=	341.04
Alignment RMSD (HBsAg)	=	0.00

Result of our Modeled Constructs

We anticipated that modifying the antibody, such as by attaching or deleting fragments, would have an impact on docking. However, we observed no changes in the docking score. This finding provides valuable insights into antibody design and modification. It suggests that the alteration of antibody structures, such as fragment attachments or deletions, may not significantly affect their binding affinity within the context of our system. This information holds particular significance as it paves the way for more flexible and versatile antibody design strategies. These findings offer new hypotheses and ideas for the design of wet lab experiments.

Mathematical Model

Mathematical models provide predictive insights into how a reaction will progress over time. With these predictions, we can adapt our experiment plans and optimize experiments by adjusting parameters such as time and concentration. This approach reduces the number of experiments needed, thereby conserving valuable time and resources.

Code availability: The code we wrote to perform the Math modelling can be found on our Software GitLab. https://gitlab.igem.org/2023/software-tools/munich/-/tree/main

In our project, we employ two types of transcytosis: Receptor-mediated transcytosis of transferrin and the melanotransferrin transcytosis pathway. Transcytosis, in general, involves several steps, including endocytosis, intracellular vesicular trafficking, and exocytosis. These processes will be described by three ordinary differential equations (ODE) and will form the foundation of our ODE model. In the model we assume antibody will not degrade.

according to figure 1 concentration change between apical chamber and the cell can be described by the concentration of antibody and the rate of endocytosis k_a and re-exocytosis k_-a concentration change between basolateral chamber and the cell can be described by the concentration of antibody and the rate of exocytosis k_b and re-endoytosis k_-b. the ODE is described as follow:

Numerical simulation

the numerical estimation is done by solving ordinary differential equation with 4th order runge-kutta method:

runge-kutta method works by estimating the slope of the solution at multiple points within a small interval and then averaging these slopes to compute a better approximation of the solution at the next time step.

h is the time step of estimation, it is set as 1 second in our simulation t_i is the initial value of time y_i is the initial value of antibody concentration k₁ is the slope calculated at the beginning of the interval. k₂ is the slope calculated at the midpoint of the interval using the k₁ k₃ is also calculated at the midpoint of the interval, using k₂ k₄ is calculated at the end of the interval, using k₃

we used estimated kinetics constants from literature[^10] ka = 1.023 * 10^-4 [s^-1] kb = 9.514 * 10^-5 [s^-1] kc = 1.358 * 10^-5 [s^-1] kd = 2.756 * 10^-5 [s^-1]

Simulation of experiment

The experimental plan of transcytosis provided by the wetlab team is as follows: First, we will apply the antibody in the apical chamber, allowing the cells to uptake the antibody. We will then wait for one hour. Afterward, all liquid will be removed, and the cell layer will be rinsed several times to eliminate any remaining antibody outside the cells. Following the rinsing process, the cell layer will be immersed in a new medium, causing the antibodies to be released from the cells. We will measure the change in concentration within the chamber over time.

import matplotlib.pyplot as plt
# Define the differential equations
def dnu_dt(nu, nc, nl, t):
    k_a = 1.023e-4
    k_b = 9.514e-5 

    dnu = - k_a * nu + k_b * nc
    return dnu

def dnc_dt(nu, nc, nl, t):
    k_a = 1.023e-4
    k_b = 9.514e-5
    k_c = 1.358e-5
    k_d = 2.756e-5  

    dnc = k_a * nu - (k_b + k_c) * nc + k_d * nl
    return dnc

def dnl_dt(nu, nc, nl, t):
    k_c = 1.358e-5 
    k_d = 2.756e-5

    dnl = k_c * nc - k_d * nl
    return dnl

# Define the 4th order Runge-Kutta method
def runge_kutta_4th_order(h, nu, nc, nl, t):
    k1_nu = h * dnu_dt(nu, nc, nl, t)
    k1_nc = h * dnc_dt(nu, nc, nl, t)
    k1_nl = h * dnl_dt(nu, nc, nl, t)

    k2_nu = h * dnu_dt(nu + 0.5 * k1_nu, nc + 0.5 * k1_nc, nl + 0.5 * k1_nl, t + 0.5 * h)
    k2_nc = h * dnc_dt(nu + 0.5 * k1_nu, nc + 0.5 * k1_nc, nl + 0.5 * k1_nl, t + 0.5 * h)
    k2_nl = h * dnl_dt(nu + 0.5 * k1_nu, nc + 0.5 * k1_nc, nl + 0.5 * k1_nl, t + 0.5 * h)

    k3_nu = h * dnu_dt(nu + 0.5 * k2_nu, nc + 0.5 * k2_nc, nl + 0.5 * k2_nl, t + 0.5 * h)
    k3_nc = h * dnc_dt(nu + 0.5 * k2_nu, nc + 0.5 * k2_nc, nl + 0.5 * k2_nl, t + 0.5 * h)
    k3_nl = h * dnl_dt(nu + 0.5 * k2_nu, nc + 0.5 * k2_nc, nl + 0.5 * k2_nl, t + 0.5 * h)

    k4_nu = h * dnu_dt(nu + k3_nu, nc + k3_nc, nl + k3_nl, t + h)
    k4_nc = h * dnc_dt(nu + k3_nu, nc + k3_nc, nl + k3_nl, t + h)
    k4_nl = h * dnl_dt(nu + k3_nu, nc + k3_nc, nl + k3_nl, t + h)

    nu_new = nu + (1/6) * (k1_nu + 2 * k2_nu + 2 * k3_nu + k4_nu)
    nc_new = nc + (1/6) * (k1_nc + 2 * k2_nc + 2 * k3_nc + k4_nc)
    nl_new = nl + (1/6) * (k1_nl + 2 * k2_nl + 2 * k3_nl + k4_nl)

    return nu_new, nc_new, nl_new

# Initial conditions
nu = 2.67  
nc = 0
nl = 0

# Time settings
t = 0
end_time = 180000
h = 1 

# Lists to store the results
times = [t]
nu_values = [nu]
nc_values = [nc]
nl_values = [nl]

# Perform the time integration
while t < end_time:
    nu, nc, nl = runge_kutta_4th_order(h, nu, nc, nl, t)
    t += h
    times.append(t)
    nu_values.append(nu)
    nc_values.append(nc)
    nl_values.append(nl)

# Print the results
for i in range(len(times)):
    print(f"Time: {times[i]} minutes - nu: {nu_values[i]} µg/mL, nc: {nc_values[i]} µg/mL, nl: {nl_values[i]} µg/mL")
# Convert the times from seconds to hours
times_hours = [t / 3600 for t in times]  # 1 hour = 3600 seconds

# Plot the results with x-axis in hours
plt.figure(figsize=(12, 6))
plt.plot(times_hours, nu_values, label='n(apical)', linewidth=2)
plt.plot(times_hours, nc_values, label='n(cell)', linewidth=2)
plt.plot(times_hours, nl_values, label='n(basolateral)', linewidth=2)
plt.xlabel('Time (hours)')
plt.ylabel('Concentration (µg/mL)')
plt.title('Concentration Changes Over Time')
plt.legend()
plt.grid(True)
plt.show()

# Find the index corresponding to 3600 seconds
time_to_find = 3600  # seconds
index_at_3600 = times.index(time_to_find)

# Get the values at 3600 seconds
nu_at_3600 = nu_values[index_at_3600]
nc_at_3600 = nc_values[index_at_3600]
nl_at_3600 = nl_values[index_at_3600]

# Print the values
print(f"At 3600 seconds - nu: {nu_at_3600} µg/mL, nc: {nc_at_3600} µg/mL, nl: {nl_at_3600} µg/mL")

Figure 2: numerical simulation of applying 2.67 µg/mL of antibody in the apical chambe, after 1 hour of incubation the concentration are n(apical):1.96 µg/mL, n(cel): 0.69 µg/mL, n(basolateral): 0.02 µg/mL respectively

Figure 3: numerical simulation with an initial concentration of 0.69µg/mL antibody in cell (obtained from figure 2), the system reached equilibrium at approximately 30 hours.

Estimation of kinetics constant with artificial neural network

The kinetic constants may differ for different transcytosis pathways and antibody constructs. To address this, we have developed a system that utilizes an artificial neural network (ANN) approach to estimate the kinetic constants[^10].

In this approach, we will collect wet lab data from transcytosis assays, which monitor changes in antibody concentration in both chambers and within the cells over time. This data will be used to train an artificial neural network consisting of 1 input layer, 1 hidden layer with 3 nodes, and 3 output nodes. The network’s input will be time, and its output will predict antibody concentrations.

Equation 3.1 describes a linear function that transforms the input time(τ) to the hidden layer. In this equation, ‘i’ represents the nodes in the hidden layer, ‘ν’ represents the weight, and ‘q’ represents the bias.

Equation 3.2 defines the sigmoid activation function used in the hidden layer, based on the values calculated in Equation 3.1.

Equation 3.3 represents the input function from the hidden layer to the output layer. Here, ‘j’ signifies the nodes in the output layer, ‘w’ stands for the weight, and ‘r’ denotes the bias.

In Equation 3.4, the sigmoid activation function is applied to the output layer using the values obtained from Equation 3.3.

Equation 3.5 defines the loss function to be optimized. In this equation, ‘p’ corresponds to the number of time points from experimental data, ‘j’ represents the number of outputs (which is 3 in this case). ‘n_j’ represents the experimental data points, such as the concentration of antibodies. In our context, ‘n₁’ corresponds to ‘n_ap’, ‘n₂’ corresponds to ‘n_cl’, and ‘n₃’ corresponds to ‘n_ba’. ‘n^ANN’ represents the predicted values from the artificial neural network (ANN).

During the training process, the ANN will work to minimize the error between the predicted antibody concentrations and the experimental data. Once the ANN is trained, we will extract the weights and biases from the neural network to calculate the kinetic constants.

Equation 3.7 represents the derivative of Equation 3.4. It serves as the sum of squared errors function to be optimized for estimating the best values for the kinetic constants. In this equation, ‘f_j’ corresponds to the three ordinary differential equations (ODEs), and ‘k’ represents the kinetic constant.By substituting Equation 3.6 into Equation 3.7, we obtain Equation 3.8. Once the optimization process is completed, the kinetic constants can be estimated. Unfortunately, due to the lack of wet lab data, this system is still in the prototype stage.

Conclusion and future prospect

Together, we used different approaches, including antibody structure prediction and docking simulations, as well as theoretical and numerical mathematical models. These methods allowed us to deepen our understanding of antibody structure, validate protein structures from two different sources through superposition, and suggest which model would likely function the best. This information can then be used to compare with wet lab results. In theory, construct W would work the best, and there is no obvious effect on docking scores from modifying the antibody, such as attaching or removing a chain or using a linker. Surprisingly, attaching an additional chain or peptide to the N-terminus does not affect the docking score, contrary to our initial belief that these additional chains attached to the variable region of the antibody might cause hindrance and affect binding. These findings provide new insights, but, of course, they also need to be further tested in the lab. The mathematical model simulation shows that the time to achieve equilibrium is approximately 30 hours, which may suggest to the wet lab team that they conduct the experiment for a longer duration.

Structure and docking simulations can be a part of antibody design. With the help of AI and other computational tools, antibody and drug development can be accelerated. Paratope-epitope prediction¹⁰ can now also be done by deep learning, and epitope-specific binding proteins have been attempted to be created by generative AI¹¹. Although there are still challenges, such as predicting structural changes upon binding, it is not surprising that complete in silico antibody design can be eventually achieved.