Click Chemistry:

The Azide + DIBAC → Triazole reaction is a powerful and versatile tool in both biological and chemical applications due to its high specificity, efficiency, and bioorthogonal nature. Biologically, the ability to selectively label biomolecules with minimal disruption is invaluable for studying complex cellular processes, developing targeted drug delivery systems, and creating stable protein-drug conjugates. Mathematically, the reaction follows a classic second-order kinetic model, where the rate of product formation depends on the concentrations of both reactants. This differential equation helps us predict how the reaction proceeds over time, providing crucial insights for optimizing reaction conditions in both research and industrial applications.

The differential equation that describes the process mentioned above is the following.

$$ \frac{d\left(1,2,3-Triazole\right)}{dt}=\left(k\cdot Azide\cdot DIBAC\right) $$

The graph describes the production of 1,2,3-Triazole using DIBAC and Azide. The y axis is in 'millimolarity’ and the x axis is in ‘time(hour)’. If we increase the parameter ‘k’ to a higher number the production will increase greatly.

Agent-Based Model

Agent-Based Model For Exosome Passage Through The Blood-Brain-Barrier:

The successful passage of the RVG29-equipped exosomes, through the BBB, is a crucial step in the drug delivery process. Due to the complex nature of the Barrier and the many interactions and changes that take place into this area, came the idea to use an Agent-Based Modelling approach in order to attempt simulating this environment, making sure that the model features key processes in a realistic as possible manner, but focusing on making a simplified yet dynamic model.

The software used for this application is Netlogo, and the reason for this choice was Agent-based modeling (ABM) , a cutting-edge computational technique, widely recognized as an exceptional tool for simulating complex natural and biological systems. It is particularly advantageous for modeling phenomena that are not readily amenable to precise mathematical formulations or traditional analytical approaches. Unlike conventional mathematical models, ABM accommodates systems characterized by heterogeneity, nonlinearity, and emergent behaviors, where strict mathematical expressions may fall short.

Elements of the ABM simulation :

Our model simulates the movement of exosomes, having taken into account the rvg-29 peptide behavior as well, within a 3D environment that represents the blood-brain barrier (BBB) and the central nervous system (CNS). The aim is to study how exosomes carrying therapeutic agents interact with the dynamic nature of the BBB. The model also tracks cytokines, including TNF-alpha, IL1-beta, and IL10, which are important for understanding inflammation responses within the system.

Setup Procedure :

The simulation begins by initializing a 3D environment, where different cell types and molecular agents are placed. The key entities in the model include:

Exosomes: These are the agents of interest in the model. They diffuse randomly in 3D space, and their movement is influenced by a combination of diffusion and the flow of blood within the artery. Also, each agent with the usage of a Normal distribution is being assigned a size value that is between 30 and 150 nanometers. Like all metrics that follow a normal distribution like the height of people, the mean is somewhere in the middle. 90 nm is exactly in the middle so that is the mean used for this Normal distribution that sets randomly the size of each exosome agent, and the values are also bounded in the said range , so they won't exceed the lower or upper bound. Lastly a probability for the possibility of exosome decay has been added, which increases at the course of time.

Endothelial Cells: These represent the cells lining the blood vessels, controlling the permeability of the BBB. They act as a barrier that exosomes must cross.

Tight Junctions: These are the structures between endothelial cells that further regulate permeability.

Pericytes and Astrocytes: These support cells of the BBB and are included in the model to represent the physical environment after exosomes pass the barrier.Astrocytes are mainly affiliated withthe health of endothelial cells, thus their permeability,while pericytes are responsible with tight junction health and neuro-inflammation regulation.

Cytokines: TNF-alpha, IL1-beta, and IL10 cytokines are introduced to model the immune response. These cytokines move randomly in 3D space and can interact with the same cellular structures. Each one has a different effect on the permeability of the barrier.

With the above being dealt with, we proceed to the next step which is the go procedure. This second and final part of the Netlogo code structure, defines the behavior of the agents, a change of their properties and how they interact with each other, wherever that is implicated.

Go Procedure :

Movement of Exosomes

Starting with the focus of the simulation, the exosome agents, they begin in the blood area, and move in an unpredictable manner. That’s to account for the diffusion of the liquid. Besides a certain change in the x-axis, each time value, due to the flow of the blood, the other two axes y and z, change almost purely due to random events. The way this movement behavior is implemented in the simulation is by using the Wiener process, for the stochastic rate change of movement for each direction.

The Wiener process is a continuous-time stochastic process, that for a set dt, it takes a value that at 99.7% of the time is between $$ (−3√𝑑𝑡,+3√𝑑𝑡) $$, this means that every set dt that passes there's a value being chosen that will belong in this range 99.7% of the time, where the range is 3 standard deviations where the standard deviation is $$ √𝑑𝑡 $$. In the simulation we have multiplied that with a constant where D=0.2(a diffusion coefficient), this is just to change D if needed and the entire term brings a decent gradual movement of the exosomes without doing big leaps in their steps. The wiener process follows a Normal distribution and has a mean equal to zero. So, the likelihood of each exosome moving to a specific direction is equal to moving to the exact opposite. This kind of behavior is also described as Brownian motion. Since we're multiplying the range with the previously mentioned constant the value each dt passes belongs in the range of $$ (−3√2𝐷𝑑𝑡,+3√2𝐷𝑑𝑡)$$.

Now that we're clear of the general movement, the behavior of the RVG-29 has also been taken into account. The way this peptide works, is that when it's located near alpha 7 acetylcholine nicotinic receptors, it tends to draw near them and proceed in attempting to bind with that receptor. Hence why the modified exosomes traverse through the barrier, with receptor mediated endocytosis. So, as each agent moves randomly in the blood area, if the random movement results in a location near the barrier, it proceeds to do small steps directed towards the barrier in order to attempt binding to the said receptor and continuing its journey inside the nervous system where the neural cells, our targets are located.

Also, the cytokines that are present in this simulation TNF-a, IL-1b and IL-10, follow this kind of random movement behavior excluding the attributes that the RVG-29 peptide adds to the exosome agents.

Interactions of agents

A)Exosomes

Now that we've described the movement of the exosomes, we can focus now on the interactions employed for this simulation. There are two interactions that the exosome agents are involved with: The endothelial cells and the neural cell targets once they have successfully passed the barrier. For a certain proximity the rvg-29 leads the exosome towards both of these cells, since both express on their surface the alpha 7 acetylcholine receptor. Regarding the interaction with the endothelial cells, since as we mentioned earlier the exosomes pass the barrier through receptor mediated endocytosis, there's a set probability that can be decreased depending on the size of the exosome that interacts. While that probability was difficult to quantify, we tried to give it a realistic value, that decreases the larger an exosome is, since it's more difficult for a larger vesicle to move past tight-junction that block large entities.

For the ones that do make it past the barrier, the continue within the CNS (central nervous system). If the random movement and leads them to a neural cell, along with the RVG-29 guidance, there's a binding probability that is very large, since the binding affinity of the RVG-29 peptide with the alpha 7 acetylcholine nicotinic receptor is well documented. That makes it highly unlikely for the exosomes not to bind. Also, the environment here is more simplified compared to that of the artery and the blood flow and many other factors that determine binding success with the same receptor on the surface of an endothelial cell.

B)Cytokines

We introduced three types of cytokines into the simulation: Tumor Necrosis Factor Alpha (TNF-α), Interleukin 1 Beta (IL-1B), and Interleukin 10 (IL-10). TNF-α and IL-1B are proinflammatory cytokines that disrupt endothelial cell function, leading to increased BBB permeability. They also downregulate key efflux proteins, such as Low-Density Lipoprotein Receptor-Related Protein 1 (LRP1) and Breast Cancer Resistance Protein (BCRP), which are responsible for clearing Aβ peptide from the CNS. This downregulation contributes to Aβ accumulation, further fueling neuroinflammation. The impact of pro-inflammatory cytokines goes beyond BBB permeability and reduced protein expression. Increased permeability allows more harmful substances, such as peripheral immune cells and pathogens, to enter the CNS, exacerbating neuroinflammation. This creates a feedback loop where higher permeability and inflammation reinforce one another.

In contrast, IL-10, an anti-inflammatory cytokine, works to counteract these effects and restore barrier integrity. However, if the levels of pro-inflammatory cytokines exceed the influence of IL-10, the cycle of increased permeability and neuroinflammation can persist unchecked.

Although we outline some of the cytokine’s roles, their actions in the simulation are simplified. Their primary effect is modeled as altering BBB permeability, with neuroinflammation building up as a result. This reflects the correlation between permeability changes and inflammatory responses, which are interdependent, with inflammation also causing a rise in permeability, thus a feedback loop is created.

C)Astrocytes and Pericytes

Astrocytes and pericytes play crucial roles in maintaining the integrity of the blood-brain barrier (BBB) and mitigating increased permeability and neuroinflammation. Positioned on the abluminal side of the barrier, these cells are essential for its proper function. Astrocytes primarily support endothelial cells (ECs) by promoting tight junction integrity, which is vital for BBB stability. In response to various stimuli, such as CNS injuries, microglial activation, and neuronal stress signals, astrocytes can produce proinflammatory cytokines like Tumor Necrosis Factor Alpha (TNF-α). This cytokine can increase BBB permeability, allowing immune cells and harmful substances to infiltrate the CNS. However, it's important to note that in the current version of our simulation, cytokine release by astrocytes has not been modeled. We are only accounting for cytokines that diffuse through the bloodstream and impact the barrier.

Pericytes, meanwhile, provide structural support to tight junctions and play a key role in regulating neuroinflammation. In response to increased permeability or neuroinflammatory triggers, astrocytes and pericytes coordinate to prevent further barrier disruption and mitigate inflammatory responses.

The Simulation

The above description outlines the functions of our simulation regarding exosome passage through the blood-brain barrier (BBB). The primary objective of developing this Agent-Based Model was to provide a clear, molecular-level understanding of the processes involved and how various parameters affect drug delivery. Each new piece of information uncovered may necessitate extensive additional hours of implementation to reflect the complexities accurately. It is important to note that while the simulation offers valuable insights, it presents a simplified scope of the procedure and serves primarily as a showcase model, illustrating fundamental concepts rather than capturing every intricate detail of the biological system.

Agent Colors in NetLogo Simulation In the NetLogo simulation, different agents are represented by distinct colors to facilitate identification: Endothelial Cells.: Red These cells are crucial for forming the inner lining of blood vessels and are represented in red. Tight Junctions: Orange Pericytes: Brown Astrocytes: Green Exosomes: Blue Cytokines (TNF-alpha): Magenta Cytokines (IL-1β): Grey Cytokines (IL-10): Yellow Neural Cells (NC): Green

The graph above illustrates the changes in endothelial cell permeability over time. As previously discussed, the modulation of permeability is influenced by factors such as astrocytes and cytokines. Additionally, pericytes play a crucial role in maintaining tight-junction integrity and reducing neuroinflammation, with the latter one contributing to decreased permeability.

Stochastic Model

A theoretical framework for the probability of an exosome engaging with the BBB

During the simulation runs of exosome passage, it became evident that exosomes must be in close proximity to the blood-brain barrier (BBB)—specifically, near the perimeter of the artery walls—to engage effectively with it. Due to the random nature of diffusion, some exosomes may never reach the barrier, which prompted the exploration of a theoretical framework for probability. This framework aims to describe the likelihood of an exosome being sufficiently close to the BBB for the RVG-29 peptide to facilitate receptor-mediated endocytosis, or at least to initiate an attempt at this critical interaction.

Regarding exosome movement in 3D space 𝑃(𝑡) = [𝑋(𝑡)𝑌(𝑡)𝑍(𝑡)] $$ X\left(t+dt\right)=\ x0\ +\ v\cdot dt\ +\ \sqrt{2D}\cdot Wx\left(dt\right) $$ $$ Y\left(t+dt\right)=\ y0\ +\sqrt{2D}\cdot Wy\left(dt\right) $$ $$ Z\left(t+dt\right)=\ z0\ +\sqrt{2D}\cdot Wz\left(dt\right) $$

Why the Rayleigh Distribution

Before diving into the probability calculations, it’s essential to understand why the Rayleigh distribution is used in the context of exosome movement. The Rayleigh distribution is ideal for modeling the magnitude of a vector in a 2D plane when its components (in this case, the height z and width y) are independent and normally distributed with the same variance, at least we assume that they have the same variance.

The Rayleigh distribution describes the distribution of the distance from the origin (in 2D space) to any point that has been sampled from two independent and normally distributed variables. In this problem, the y-axis (width) and the z-axis (height) represent such independent, normally distributed random variables because the exosome’s movement in the artery (in these dimensions) is stochastic and driven by diffusion. In our model, although the exosome moves in 3D space, the longitudinal movement (along the x-axis) is not relevant for determining the likelihood of engaging with the BBB. Thus, we focus on the 2D plane formed by the y and z axes. The Rayleigh distribution fits this scenario as it allows us to calculate the Euclidean distance (or radius) in 2D space, which is the key factor determining whether the exosome is close enough to the barrier. The probability that the exosome is beyond a certain distance (or radius) from the BBB can be calculated using the cumulative distribution function (CDF) of the Rayleigh distribution. This helps describe the likelihood that the exosome reaches the desired proximity to engage with the BBB.

So, the radius of the exosome at any moment is equal to $$ r\left(t+dt\right)=\sqrt{z\left(t+dt\right)^2+y\left(t+dt\right)^2} $$

We aim to calculate the probability that, at any given moment, the radius R exceeds a threshold r0, where r0 is a lower bound guaranteeing that an exosome can engage with the BBB. This gives the following probability expression: So, we want the probability that $$ R>r0-r\left(t\right)\rightarrow P\left(R>r0-r\left(t\right)\right)=1-P\left(R\le r0-r\left(t\right)\right) $$

$$ P\left(R>r0\right)=1-\exp \left(-\frac{r0^2}{2σ^2}\right) $$ with $$ σ=t $$

However, this standard CDF does not account for the continuous change in the exosome’s position due to the random nature of diffusion. To incorporate the ongoing updates to the exosome’s location over time, we need to modify this CDF using the Wiener process. At each small-time increment dt, the Wiener process updates the values of y and z, which means we must adjust the variance based on the sampling frequency that updates the position of the exosome. The variance of each Wiener process at time dt is 2D⋅dt, where D is the diffusion coefficient. Thus, the time-dependent variance for the Rayleigh CDF becomes:

$$ \sigma =\sqrt{2\cdot D\cdot dt} $$. So $$ 2\cdot \sigma ^2=4\cdot D\cdot dt $$

To describe the evolving probability of an exosome being within the desired radius r0, we now modify the Rayleigh CDF to reflect the dynamic position of the exosome:

$$ P\left(R>r0-r\left(t+k\cdot dt\right)\right)=1-\exp \left(\ -\frac{\left(r0-r\left(t+k\cdot dt\right)^2\right)}{2\sigma ^2}\right) $$

And with $$ 2σ^2\ =4\cdot D\cdot dt $$

$$P\left(R>r0-r\left(t+k\cdot dt\right)\right)=1\ -\ \exp \ \left(\ -\frac{\left(r0-r\left(t+k\cdot dt\right)\right)^2}{4\cdot D\cdot dt}\right) $$

Where $$k\ ∈\ ℕ $$ and $$ \ ℕ (1,2,3,...) $$

Molecular Docking

Molecular Docking and Dynamics of RVG-29 and α7-nAChR

One of the key elements in enhancing the drug delivery of exosomes carrying microRNA195 is the ability of the RVG-29 peptide (Rabies Virus Glycoprotein) to specifically target the α7 nicotinic acetylcholine receptor (α7-nAChR). This interaction improves the efficiency of exosome-mediated delivery into the central nervous system. To explore this binding interaction, we conducted molecular docking and molecular dynamics simulations to confirm the stability and strength of the RVG-29 and α7-nAChR complex. For the receptor structure, we obtained the crystal structure of the ligand-binding domain of α7-nAChR from the Protein Data Bank (PDB). The original PDB file contained two α7-nAChR receptors, but we removed one to focus exclusively on a single receptorligand interaction. Since the RVG-29 peptide was not available in the PDB database, we turned to the AlphaFold server by Google DeepMind. AlphaFold uses advanced AI to predict the 3D structure of proteins with high accuracy based on their amino acid sequences. This allowed us to generate the structure of RVG-29 for use in the docking simulations.

Docking Simulation

We employed HDOCK, a web-based platform specialized for protein-protein docking, to simulate the binding of RVG-29 to α7-nAChR. The key results from the docking simulation were as follows: • Docking Score: -300.94 • Ligand RMSD: 42.63 • Confidence Score: 0.8164 While the docking score indicates a strong interaction (lower scores suggest stronger binding), the RMSD value is relatively high. RMSD measures the structural deviation between the docked pose and the reference structure, with lower values generally indicating better stability and accuracy. However, it is important to note that RVG-29 is a flexible peptide with multiple rotatable bonds, and docking algorithms sometimes struggle to capture this flexibility, which can lead to artificially high RMSD values. In cases like this, RMSD alone may not be the best measure of binding quality.

Molecular Dynamics

Molecular Dynamics Simulation

To further assess the stability of the RVG-29 and α7-nAChR complex, we performed Molecular Dynamics (MD) simulations using GROMACS, a widely-used and powerful MD simulation tool. Given the computational intensity of MD simulations, we utilized the resources of the High-Performance Computing Laboratory (HPC Lab) at CEID, University of Patras. MD simulations provide a more realistic view of molecular interactions by allowing both the ligand and receptor to move and adjust naturally over time. This helps to correct any inaccuracies in the initial docking pose and provides insights into the long-term stability of the complex. During the simulation, we observed that the RVG-29 peptide remained stably bound to the ligand-binding domain of α7-nAChR, supporting the reliability of the initial docking results despite the high RMSD.

Force Field Selection: The force field chosen for the molecular dynamics between the two chosen proteins is the CHARMM36m chosen for protein-to-protein docking and the tip3 water model was used as the solvent for the simulation.

Simulation Box: A dodecahedral simulation box was constructed around the molecule of interest, positioning it at the geometric center of the box. The box was designed to ensure a minimum distance of 1.0 nanometer (nm) between the molecule and the boundaries of the box, providing sufficient space for solvation and reducing interactions with periodic images. The use of a dodecahedral shape was chosen for its computational efficiency, as it requires less solvent than a traditional cubic box while still maintaining periodic boundary conditions. To neutralize the system's overall charge, appropriate numbers of sodium (Na⁺) and chloride (Cl⁻) ions were added to the simulation, ensuring proper electrostatic balance and enabling realistic biological simulations.

Energy Minimization: Energy minimization was carried out using the steepest descent algorithm, with a maximum force tolerance (emtol) of 1000.0 kJ/mol/nm. This step aimed to relieve steric clashes and stabilize the system by minimizing unfavorable interactions. A step size of 0.01 nm was applied, and the minimization was allowed to run for up to 50,000 steps to ensure convergence. For interaction calculations, the Verlet cutoff scheme was used for neighbor list updates, while electrostatic interactions were treated with the Particle Mesh Ewald (PME) method, applying a 1.0 nm cutoff for both electrostatic and Van der Waals interactions. Periodic Boundary Conditions (PBC) were imposed in all directions, simulating an infinite system. This energy minimization step ensured the system was at a low-energy state, free of major conflicts, and ready for molecular dynamics simulations.

Equilibration: The equilibration process was performed in two phases: NVT (constant Number of particles, Volume, and Temperature) and NPT (constant Number of particles, Pressure, and Temperature).

NVT Equilibration: The NVT simulation aimed to stabilize the temperature of the system. Using a leap-frog integrator, we conducted 50,000 steps (100 ps) with a 2 fs time step. The protein was position-restrained with -DPOSRES to maintain structural stability. The V-rescale thermostat regulated the temperature at 309.65 K, with a coupling time constant of 0.1 ps. Initial velocities were assigned from a Maxwell distribution at the target temperature, and 3D periodic boundary conditions were applied to prevent artifacts from finite size.

NPT Equilibration: Subsequently, we performed NPT equilibration for another 50,000 steps at the same time step. The Parrinello-Rahman method managed pressure coupling, maintaining a reference pressure of 1.0 bar with a 2.0 ps time constant. The protein remained positionrestrained to avoid distortion. Long-range electrostatics were calculated using Particle Mesh Ewald (PME) with a cutoff of 1.0 nm for both Coulomb and Van der Waals interactions. This two-step equilibration ensured the system achieved the desired temperature and pressure, providing a stable foundation for subsequent production simulations.

In the visualization produced by VMD, the RVG-29 peptide is prominently displayed at the center as a multi-colored structure. The similarity in color between certain parts of the receptor and the peptide arises from the software's selection process. VMD allows us to choose a specific sequence of residues, which results in the selection of the first 29 residues across all chains. This includes the five chains of the alpha-7 nAChR, leading to overlapping color representations.

The hydrogen bond graph derived from the GROMACS simulation over 10,000 picoseconds indicates that the system is really stable, beginning with over 49,200 hydrogen bonds at the initial docked state. Notably, there is a sudden increase in the number of hydrogen bonds at the start of the simulation, reaching a peak of approximately 49,600 bonds. This strong bonding affinity is consistently maintained throughout the simulation duration. While the graph reflects the total hydrogen bonds within the system, the initial spike can be attributed to the high binding affinity between RVG-29 and the alpha-7 nAChR. Consequently, the fluctuations in the total number of hydrogen bonds are influenced significantly by the interactions between these two proteins. Overall, the dynamics observed in the simulation affirm a robust bond and high attraction between RVG-29 and the alpha-7 nAChR.

Future Steps

Beyond the iGEM competition, we are eager to continue refining our Agent-Based Model to explore various factors that connect Alzheimer's disease with the functions of the Blood-Brain Barrier (BBB). One promising avenue for future research is the relationship between the TNF-alpha cytokine and its effects on the BBB, as it may play a critical role in the progression of the disease. Agent-Based Models can provide valuable insights into the complex interactions and mechanisms underlying diseases like Alzheimer's. The knowledge generated from these models could be instrumental in improving our understanding of the disease and potentially informing new therapeutic strategies.