Cycle 1: Design

Since one of the parts of our project is to create a fusion protein system for the bottleneck reaction of naringenin synthesis, we decided to model it in silico. The model’s primary purpose is to check whether the shorter distance between active sites leads to the higher production of naringenin. In the process, we found out that even the modeling part required to cover the main steps of engineering: design, build, test, and learn. Therefore, we began with the already existing modeling workflow described for a fusion construct with a malaria pre-erythrocytic vaccine candidate [1]. The described modeling approach was initially designed as follows:

1. Choosing the protein candidates
2. Selecting linkers
3. Primary structural analysis of the fusion system
4. Protein modeling
5. Scoring modeled complexes to choose the model for visualization
6. Plotting Ramachandran graphs before structure refinement
7. Structure refinement
8. Plotting Ramachandran graphs after structure refinement

Cycle 1: Build

The building stage of the modeling flow consisted of performing the initial design flow. In this cycle we have raised a hypothesis that homology-based protein modeling approaches should not be the most suitable method for our case, since there are few fusion models existing.

1. The protein candidates were chosen to be 4-coumarate-CoA ligase 2 from Glycine max and chalcone synthase from Arabidopsis thaliana
2. Seven linkers were chosen for our fusion system. Four of them were flexible glycine and serine linkers (GSG, (GGGGS)n, where n is equal to 1, 2, and 3) and the other three were rigid ((EAAAK)n, where n is equal to 1, 2, and 3)
3. For a primary structural analysis we studied PDB: 3TSY structure - an experimentally determined structure of the fusion system that has 4CL protein in it. We took this structure as a starting point of what we can expect from our models
4. Protein modeling was initially attempted to perform using homology-based modeling program SWISS-MODEL[2]
5. The modeled structures in the first cycle were evaluated using VoroMQA[3], QMEAN[4], and QMEANDisCo[5] scores
6. The plot was not drawn in this cycle
7. The structure refinement was not performed in this cycle
8. The plot was not drawn in this cycle

Cycle 1: Test

The output of SWISS-MODEL showed that sequence identity between the templates found and the sequence that is modeled is 66.42 - 69.36%. According to the specialists in protein modeling, homology-based methods are suitable, when the sequence identity between the template and the modeled sequence is more than 25%, belonging to the “daylight” zone[7], thus, in theory, the method could be applied.

The modeling flow was tested with PyMOL[6] visualization program using commands `cealign`, `set seq_view`. We were looking for an accurate representation of the linked proteins in the fusion system by comparing them to their distinct versions. The main focus on the system was the linker region, which in the homology-based modeling case was composed of the uncoiled domain of 4CL that made up a highly disordered massive linker region. We stopped the flow after the visualization because we decided that structure refinement using molecular dynamics will not introduce any significant changes to the structure.

Cycle 1: Learn

Nonetheless in theory the homology-based modeling method could be applied, these modeling results proved that this approach could not be applied to our fusion protein system. Therefore, we decided to try running ab initio protein modeling programs.

Cycle 2: Design

The design (main steps) of the modeling flow stayed the same as it was in the first cycle, yet also we included one more step - energy minimization with Yasara[10] - that was inserted after the protein modeling step. The main changes in the second cycle were the choice of protein modeling programs. In this engineering cycle we took ab initio modeling programs trRosetta[8] and RoseTTAFold[9] to model our proteins.

Cycle 2: Build

1. Protein candidates stayed the same as in the first cycle
2. Linkers stayed the same as in the first cycle
3. Primary structure analysis was not required to be redone in the second cycle
4. Protein modeling was performed using trRosetta and RoseTTAFold simultaneously
5. The modeled structures were evaluated using VoroMQA score
6. Energy minimization with Yasara.
7. The Ramachandran plots were drawn before molecular dynamics (MD) simulations
8. We performed short MD simulations for the trivial linker (GSG) case
9. The Ramachandran plots were drawn before molecular dynamics (MD) simulations

Cycle 2: Test

The modeled structures visually were not satisfactory - the disordered domain of 4CL protein was still present in the structure. The energy minimization step with Yasara and MD simulations did not introduce any significant changes to the structure.

Cycle 2: Learn

After this step we decided to consult the bioinformatician who works with protein modeling. He suggested we include multiple sequence alignment files into our modeling flow. Additionally, we got advice to use more scoring functions in the evaluation step. Yet also we consulted the protein molecular dynamics specialist to get a professional insight into how we should run the MD for our system to get a more significant output. After the consultation we adjusted the box size of the system and the duration of the simulations.

Cycle 3: Design

In the third engineering cycle we included our own generated multiple sequence alignment (MSA) files as input to the protein modeling programs of our choice. In addition, we took into account the advice from the specialists and included QMEAN and QMEANDisCo scores into the evaluation of the structures. The structure refinement was not performed in this cycle due to the lack of computational resources at the time.

Cycle 3: Build

1. Protein candidates stayed the same as in the previous cycle
2. Linkers stayed the same as in the previous cycle
3. Primary structure analysis was not required to be redone
4. Protein modeling was performed using RoseTTAFold
5. The modeled structures were evaluated using QMEAN, QMEANDisCo, and VoroMQA scores
6. The plot was not drawn in this cycle
7. The structure refinement was not performed in this cycle
8. The plot was not drawn in this cycle
9. The distances between active sites in fusion protein systems with rigid linkers were calculated

Cycle 3: Test

The modeled structures were visually satisfactory - the disordered domain of 4CL protein was less disordered in the modeled structure.

Cycle 3: Learn

The multiple sequence alignment (MSA) files have a significant impact on the modeling output.

Cycle 4: Design

According to the output of the third cycle, the fourth engineering cycle was not required. However, after the latter cycle was finished, AlphaFold2[11] became available for public use. Therefore, we decided to apply this highly evaluated tool in our modeling workflow. The structure refinement was not performed in this cycle due to the insignificant impact of MD for structures recorded in the studies[12].

Cycle 4: Build

1. Protein candidates stayed the same as in the previous cycle
2. Linkers stayed the same as in the previous cycle
3. Primary structure analysis was not required to be redone
4. Protein modeling was performed using AlphaFold2
5. The modeled structures were evaluated using QMEAN, QMEANDisCo, ProQ2D, ProQRosCenD, ProQRosFAD, ProQ3D scores
6. The plots were drawn in this cycle
7. The structure refinement was not performed in this cycle
8. The distances between active sites in fusion protein systems with rigid linkers were calculated

Cycle 4: Test

The output of the fourth cycle had the less disordered domain just as the output of the third cycle.

Cycle 4: Learn

We were visually satisfied with the modeled complex. Nonetheless, AlphaFold2 models do not require the usage of molecular dynamics for structure refinement, we considered applying them to get a better insight of the system with flexible linkers.