key: cord-0846939-w3c5s66o
authors: Xu, Zhaobin; Zeng, Qiangcheng
title: More or less deadly? A mathematical model that predicts SARS-CoV-2 evolutionary direction
date: 2022-03-10
journal: bioRxiv
DOI: 10.1101/2022.03.10.483726
sha: 3a13db58d27ba2cc3358ab3e90fbd3109f8571c8
doc_id: 846939
cord_uid: w3c5s66o

SARS-CoV-2 has caused tremendous deaths world wild. It is of great value to predict the evolutionary direction of SARS-CoV-2. In this paper, we proposed a novel mathematical model that could predict the evolutionary trend of SARS-CoV-2. We focus on the mutational effects on viral assembly capacity. A robust coarse-grained mathematical model is constructed to simulate the virus dynamics in the host body. Both virulence and transmissibility can be quantified in this model. The relationship between virulence and transmissibility can be simulated. A delicate equilibrium point that optimizing the transmissibility can be numerically obtained. Based on this model, we predict the virulence of SARS-CoV-2 might further decrease, accompanied by an enhancement of transmissibility. However, this trend is not continuous; its virulence will not disappear but remains at a relatively stable range. We can also explain the cross-species transmission phenomenon of certain RNA virus based on this model. A small-scale model which simulates the virus packing process is also proposed. It can be explained why a small number of mutations would lead to a significant divergence in clinical performance, both in the overall particle formation quantity and virulence. This research provides a mathematical attempt to elucidate the evolutionary driving force in RNA virus evolution.

1 Introduction： 26 History indicates that no lethal virus pandemic last forever. For example, the 1918 flu, which claimed tens 27 of millions of lives, was evolved into a significantly less deadly seasonal flu after 1920 [1-2]. Similar 28 situation happened such as 2009 H1N1 virus gradually lose its virulence soon after the pandemic [3] [4] . The 29 pandemic caused by RNA virus is different from those caused by bacterium in a way that it normally 30 engages a quick virulence declination after the outbreak. As far as scientist can tell, the bacterium that 31 caused the black death never lost its virulence [5] . Now the question emerges: will SARS-CoV-2, the virus 32 that causes COVID-19, follow a similar trajectory? So far as we observe, there is a significantly declination 33 of SARS-CoV-2's virulence. For instance, the virulence of omicron strain is less than the delta strain, and 34 the virulence of delta strain is less than the original strain. However, lots of people concerns the emergence 35 of a more virulence strain. It is useful when we refer to the history but history itself is also highly unreliable. 36 Therefore, it is necessary to derive a scientific theory that can help predict the evolutionary direction of 37 SARS-CoV-2. 38 The scientific exploration toward this question started in early 1980s, when two mathematical biologist Roy 39 Anderson and Robert May, proposed that germs transmit best when hosts shed a lot of the pathogen [6-7]. 40 virulence and transmissibility go hand in hand, until the germ gets so deadly it winds up killing its host too 41 soon, and therefore can't spread at all. This is known as the transmission-virulence trade-off [8] [9] [10] . The own strategies to spread, and some of those strategies allow the germ to maintain high 53 virulence and transmissibility. 54 Here we proposed a new mathematical theory that predicts the virus evolutionary direction. Instead of using 55 logic language and examples, we developed a mathematical approach that predict the evolutionary trend of 56 RNA virus. We focus the effect of mutation on viral assembly capacity and explain why the alternation of 57 viral assembly capacity will significantly change its virulence and transmissibility. 58 The viral assembly capacity has been extensively studied [13] [14] [15] [16] [17] [18] [19] [20] . However, little attention is paid to linking 59 this assembly kinetics into its virulence evolution. In the result section 3.1, we first illustrate why and how 60 this self-assembly capacity influence its replication and the formation of virus particles. In 3.2 we apply this 72 We used ordinary differential equations to describe the dynamics of the virus together with antibodies. 14 73 reactions are constructed: Reaction (1) is the replication process of virus, reaction (2) and (3) are the translation process of virus 124 mRNA, reaction (4) is the packaging process of virus, reaction (5)-(8) is the interaction process between 125 antibody and virus, and reaction (9)-(12) is the degradation process of different components of the virus.

Reaction (13) to (14) are the antibody regeneration processes. We think that the evolution track of the virus 127 is the process of maximizing the production of virus particles because only complete virus particles are 128 infectious, and neither single mRNA nor protein is infectious. The ordinary differential equations based on 129 those reactions are listed below: 

( 1) = − 1 + − 1 ;

( 2) = − 2 + − 2 ;

(25)

Equation (24) was not taken into consideration. The proliferation of the antibody was linearly correlated to its binding 133 complex since this binding complex will further stimulate the regeneration of specific antibodies. We treat 134 the immune clearance of antibody binding complex is very fast and efficient. Therefore, k7 and k9 was 135 assigned to 1 in the following simulation process. Meanwhile, we proposed that the degradation rate was 136 significantly different between mRNA and proteins. We assigned a larger value to k10 while keeping k11 as 137 a small number. To better illustrate how mutation can dramatically influence its toxicological characteristics, a kinetic model 141 that simulating the virus packing process is constructed. The agglutination nucleus is treated as the mRNA.

In general, the virus packing process can be represented as multiple capsid proteins assembling around their 143 core mRNA. This would eventually lead to the formation of a complete virus particle. The virus assembly 144 process is a highly complicated process which might include the induction of signaling molecules, the 

.

An +B AnB (forward constant = c1; reverse reaction constant = 3 1 −1 ) 176 (26) 177 178

The ordinary differential equations based on those reactions are defined as below: 179 180 ( ) = −2 * 1 * * + 2 * 2 * 2 − 1 * ( 2 + 3 + ⋯ . . + −1 ) * + 2 * ( 3 + 4 + ⋯ + ) − 1 * ( + + 2 + 181 ⋯ + −1 ) * + 3 * + 2 * 3 1 * ( 2 + ⋯ + );

(27) 182

( 2 ) = 2 * 1 * * − 2 * 2 * 2 − 1 * 2 * + 2 * 3 + 3 2 1 * 2 − 1 * 2 * ;

(28) 184 185 ( ) = 1 * −1 * − 2 * − 1 * * + 2 * +1 − 1 * * + represents the final packed virus particle which has n structure proteins and one 202 mRNA inside. represents the incomplete packing intermediates which has i structure proteins and one result section 3.1, we will explain in detail how we determine the reaction coefficient of each reaction in the 205 process (26). also weakly affected by mutations. Therefore, reaction (4) is the most sensitive process toward mutations.

This will be fully illustrated in our second model. Mutations will lead to a slight increment in the binding 231 energy between capsid protein and mRNA, they might also lead to a slight increase among capsid monomers.

This slight change in binding kinetics will lead to a significant difference in the final formation of complete 233 virus particle. The distinguished point of our theory is that we proposed that the enhancement of packaging 234 ability will greatly weaken its replication ability, because we can see from reaction (1) that only the naked 235 mRNA can make further replication, while the mRNA wrapped by capsid protein cannot make further 236 replication and translation. Therefore, the improvement of virus packaging ability will inevitably lead to the 237 decline of virus replication and translation ability. It is the increase in the number of virus particles that leads 238 to the decline of its reproductive capacity, followed by the decline of the overall mRNA level and the 239 content of all translation proteins, which also corresponds to a relatively low immune response (lower 240 antibody production). We believe that the toxicity of virus depends on the total amount of heterologous 241 substances produced, that is, the sum of mRNA, replicase and structural protein. This toxicity can also be 242 measured by the highest level of antibody generated by stimulation. With the improvement of virus 243 assembly ability, the number of virus particles will be significantly increased, and the infectivity will be 244 significantly enhanced, but its toxicity will be significantly reduced. However, this assembly ability will not 245 be enhanced indefinitely. Take two extreme cases for example. The first extreme case is that the virus 246 packing capacity is 0, which means that it is impossible to package into virus particles. In this case, the 247 number of virus particles is naturally very small, almost zero. In other extreme condition when the assembly 248 ability is super strong, and all mRNA will be wrapped by its structural protein to form virus particles. In this case, because there is no naked mRNA, the virus cannot be further copied and translated, so the number of 250 virus particles finally formed is also very small, almost zero. Therefore, there is an optimal binding force or 251 binding coefficient k5 and k-5 between these two extreme cases, which can ensure the maximum virus 252 particles produced. We think that this is the driving force behind the evolution of virus. The original 253 COVID-19 strain has a weak packaging ability, so it replicates faster and is highly toxic. Meanwhile, its 254 virus particles are in low quantity, thus has a relatively low transmissibility comparted to late strains. With 255 the extension of time, the strain with strong packaging ability will gradually replace the strain with weak 256 packaging ability, because the number of complete virus particles produced by it will increase, but this 257 packaging ability will not increase indefinitely, and there is an optimal packaging ability to produce the 258 largest number of virus particles. In result section 3.2, we will illustrate how the change in packing capacity 259 will influence the virus evolution.

To better demonstrate why tiny binding energy change between capsid monomer and mRNA, together with 261 the slight improvement in binding energy among those capsid monomers, would lead to a significant change 262 in the kinetic coefficient in reaction (4) should be c2/c1*c3*c1. Since we treat the forward reaction constant as c1, the reverse reaction kinetic 287 coefficient can be derived to be c2*c3/c1.

For the reaction 290

Ai +B AiB, the binding energy is i times of the binding energy between A and B. Therefore, the 1 . Since we treat the forward reaction kinetic as c1, the 292 reverse reaction kinetic coefficient can be derived as 3 1 −1 . The assembly process of SARS-CoV-2 is 293 normally contributed by hundreds of capsid proteins. Therefore, n should be a very large number. Here for 294 simplicity, we use n = 5, 10 and 20 to perform the simulation in section 3.3. We illustrate the bigger the n be, 295 the more significant the scale effect would be. A tiny alternation of binding energy among monomers could 296 lead to a dramatic change in the final particle formed. k1 = 1e-5; k-1 = 1e-5; k2 = 0.15; k3 = 0.5; k4 = 0.5; k5 = 1e-3; k-5 = 0.9; k6 = 5e-6; k-6 = 1e-14; k7 = 1; k8 = 5e-6; k-8 = 1e-14; k9 = 1; k10 = 0.1; k11 = 0.002.

We simulate the virus particle dynamics together with antibody dynamics in figure 1. As shown in figure 1 , 319 when the mutation increases the viral assembly capacity, the reverse reaction constant k-5 will decrease. The 320 decrease in k-5 will lead to a smaller quantity of antibodies, which corresponds to a milder virulence. The 321 overall particle number can be calculated as the cumulative virus number over time. Although the peak virus 322 particle loading amount of the dashed yellow line is smaller than that of the dashed blue line, its cumulative 323 number over time is larger than that of the dashed blue line. Therefore, mutations that promote the virus 324 packing process will display an evolutionary advantage over its original strain. However, this trend is not 325 continuous. As shown in the dashed green line, when the virus packing capacity surpasses certain threshold, 326 a further decrease in k-5 will lead to a smaller cumulative number. The logic behind the trade-off between 327 virulence and transmissibility in our model is that we assume the virus would optimize its shedding particle 328 scale. Any improvement in the packing capacity will form an obstacle to its RNA replication and translation, 329 thus decrease the overall heterogenous proteins. This will display a reduction in virulence in clinical 330 performance. Before reaching the optimal point of maximizing the overall particle number, the improvement 331 in packing capacity will lead to the generation of a more completed packed virus. This will display as a 332 strong enhancement in transmissibility with a larger reproduction R0 value. After surpassing this optimal 333 point, further mutations lead to faster packing will lose their evolution advantages since fewer virus particles 334 would be generated. Therefore, the declination of virulence is not perdurable. The virulence of specific virus 335 might engage significant declination at the early stage of the epidemic but would remain at a certain range 336 after it reached the delicate balanced point. We proposed that this might be the underlying mechanism 337 behind the intricate trade-off relationship between virulence and transmissibility. k5 can be calculated as 1e-6*1.05^(i), i is the value displayed in x-axis. k-5 can be calculated as 1*0.95^(j), j 346 is the value displayed in y-axis. 347 It can be seen from figure 2, at the same k5 value, the overall virus particles generated are very sensitive 348 toward the dissociation constant k-5. Mutations that promote viral assembly will lead to a significant 349 increase in the overall virus formation quantity. When j increases in y-axis, k-5 will decrease 350 correspondingly, the overall virus quantity will engage a significant growth before it reaches an optimal 351 point. The forward reaction constant k5 also influences the overall virus particles number, but not as 352 significant as k-5. The binding energy can be transformed into the equilibrium constant kd, which is equal to 353 k-5/k5. We could obtain an optimal combination of k5 and k-5 using a global searching approach. We could 354 also calculate the optimal kd value which is directly correlated with binding energy. We can see from figure 3a, the optimal value of k-5 will increase under a higher initial antibody level. This 367 means massive vaccinations and infections will promote the virus to evolve into a more toxic strain. The 3.2.4 Antibody binding features in different species will cause different optimal points in different species.

Besides the antibody levels studies in the previous section, the antibody attribute would also significantly 379 impact the final optimal equilibrium point. We supposed the binding capacity has a positive relationship 380 with k-5 value. In species with fast binding kinetics, their k6 and k8 value are more prominent than human 381 beings. In this case, the evolution process will select the strain with weak assembly capacity with an 382 enormous k-5 value. The assembly capacity tends to be strongly truncated in species with suitable antibody 383 binding kinetics such as bats. We gradually recognized the unique aspects of antiviral immune response in . Therefore, bats should have a bigger k6 and k8 value than human beings. We also evaluate the 389 influence of k-6 and k-8 on the optimal packing point, and we found the consequences of those reverse 390 dissociation constants are very tiny. The modeling results are different from our expectations. As shown in figure 4a, the decrease of k6 could disfavor the virus packing process, enhancing its dissociation rate. The 392 optimal dissociation constant k-5 will increase first, as shown in the green cycle marked as region 1. After it 393 surpasses a threshold, it would display an opposite trend, as shown in the red area marked as region 2. We 394 predict the cross-species transmission of coronavirus such as SARS, MERS, and SARS-COV-2 in region 2.

The transmission from strong immunity species to weak immunity species will lead to a fatal pandemic in The coarse-grained model is just a simplified representation of an actual situation. Although it reflects the 425 subtle relationship between virulence and transmissibility, people might argue that it is not reasonable to 426 have such dramatic change in k-5 value based on thermodynamic principles. For instance, according to the 427 equation ln kd = -ΔG/(RT), the ΔG only engages minor changes due to mutation. Therefore, kd value and 428 k-5 value cannot be altered in large magnitude. We have to reiterate that the coarse-grained model is not 429 physically linked to the kinetics behaviors. We developed a small-scale model that can be physically 430 connected to the binding energy perturbation due to mutation to simulate the viral assembly process better. The self-assembly of a set of viral proteins into a regularly shaped capsid involves a fine-tuned range of 433 interactions between the capsid proteins to create this highly ordered, symmetrical, supramolecular structure.

In a simplified view, assembly is a spontaneous process driven by weak protein-protein interactions on the CPs, which arrange as trimeric subunits. The whole range of particle sizes was found in this study [38] , from 454 complete capsids, via large incomplete capsids, to small capsid intermediates ( Figure 5 ). Full capsids were 455 analyzed separately according to their capability to be permeable to uranyl ions, as characterized by TEM 456 experiments. The right graph in Figure 5 shows the assembly evolution of intermediates to complete capsids 

Since the experiment in figure 5 only uses capsid proteins, we also removed mRNA in our modeling. As 469 shown in figure 5, the complete packed particle will gradually increase through time. This trend can also be 470 reflected in our model. The completely packed particle is represented as greed line in figure 6 . It can be seen 471 that the fully packed virus particle will also increase through time. Beside that, the poorly packed particles After the validation process, we try to explain why a few mutations would cause a significant alternation in 480 the virus formed. We will illustrate how the number of capsid proteins in one particle influences the overall 481 packing capacity. Three cases were compared in our study when we utilized n = 5,10,20 respectively. In all 482 of those three cases, we think mutations would only lead to a minor enhancement toward the binding energy. 

It can be seen from Figure7A that the mutation only causes 10% enhancement in binding affinity. The 496 mutant strain displays a more robust particle formation capacity with a 50% increase than the original strain.

When n increases, the magnification would be more significant as shown in figure 7b and figure 7c . For 498 example, this enhancement was more than 100% when n = 10; and was 400% when n = 20. As we described 

The declination in virulence might also be explained using the second model. It was noticed the newly 509 emerged strain, such as omicron, displays weak virulence. It is hard to explain why small mutations would 510 lead to a significant decline in virulence, especially in the early epidemic phase from March 2020 to May 511 2020. The overall death rate experienced a significant decrease, from 10% to 1% [41] . Although we cannot 512 deny the possibility from other aspects such as UTR deletion proposed in our previous preprint publication 513 [42], we explain that point mutation might also contribute to the rapid declination in virulence. The 

It can be seen from figure 8 that the concentration of naked mRNA is decreased by 19% after mutation.

Although this 19% is not significant, those naked mRNA will participate in further replication and 521 translation. Therefore, a noticeable concentration difference in the equilibrium state might lead to a 522 substantial difference after several rounds of replication. In the discussion part, we discuss the possible application of this theory, especially toward those sensitive 527 questions listed in the introduction part. 528 Ⅰ:What is the exact relationship between virulence and transmissibility? 529 According to our coarse-grained model, there would be a balance point between virulence and 530 transmissibility. The evolution will advance the formation of mutants that generate more virus particles. The 531 key feature of our coarse-grained model is that we focus on the virus packing process on the overall virus 532 formation. It is rarely noticed before that the strong assembly kinetic might inhibit the proliferation of new 533 mRNA and proteins. The physical background of this model is that we think the mRNA is only active in the 534 non-binding state. The binding of the capsid will hamper its replication and translation. Our simulation 535 results indicate that the mutant with increased assembly capacity will be selected given an impoverished 536 packing background. The virus shedding amount will engage significant enhancement, and its 537 transmissibility will boom, as displayed in the case of SARS-CoV-2. 538 Meanwhile, as more mRNA are trapped around capsid protein, its replication capacity will decrease, 539 displaying a weak virulence in the clinic. This trend can also occur in other RNA virus infections such as the 540 process, is not a historical coincidence. This question is critical since it greatly influences our prevention policy. While some scientists were 545 concerned that a fatal mutant would eventually evolve [43] , our model rejects this possibility. There might 546 be some fatal mutants emerging, but the evolution pressure will limit the large-scale transmission of those 547 toxic strains. Therefore, the overall evolutionary trend of SARS-CoV-2 is to deduce its virulence. Only by 548 reducing their virulence can they increase their transmissibility. Consequently, we should not worry too 549 much about the emergence of a more deadly strain since it is strongly impeded in evolution. However, as we In section 3.3, it is mathematically demonstrated that a few mutations can lead to a significant change in the 574 overall particle formation together with its virulence. While many researchers focus on the effect of 575 mutation on spike protein which directly influences the virus entrance capacity, we suggested modifications 576 on capsid protein or mRNA might play an essential role in determining its clinical features. Mutations that 577 promote the binding among capsid proteins or the binding between capsid and mRNA would significantly 578 enhance the final formation of complete particles. The polymeric aggregation process will amplify the tiny 579 binding energy change into a tremendous difference in its assembly outcome. Therefore, it does not 580 necessarily need many mutations to change its clinical severity dramatically. We should not evaluate its 581 toxicity based on its homology and similarity to the original strain. For example, very high homology exists In the end, we have to admit there are many uncertainties in our model. This model is just a simple 585 mathematical representation of a real scenario. The parameter values can be further improved and modified 586 given sufficient experimental data. It needs further experimental or ecology validation eventually. 587 Meanwhile, we did not intend to encourage people to reduce protection against SARS-CoV-2 when we 588 acclaim the virulence of SARS-CoV-2 will engage a declination. This work only provides a hypothesis on 589 the evolutionary relationship between virulence and transmissibility. 

Federal Reserve Bank of 592 New York, 2020, 5

Patterns of mortality during pandemic: An example of Spanish flu pandemic of 1918

Overview/reflections on the 2009 H1N1 pandemic

Pandemics, places, and populations: Evidence from the Black Death

Population biology of infectious diseases: Part I

Coevolution of hosts and parasites

Transmission-virulence trade-offs in vector-borne diseases

Virulence evolution and the trade-off hypothesis: history, current state of affairs and the 604 future

The virulence-transmission trade-off in vector-borne plant viruses: a review of 606 (non-) existing studies

Next step in the ongoing arms race between myxoma virus and wild rabbits in Australia is 608 a novel disease phenotype

Pathogen survival in the external environment and the evolution of virulence

Architecture of African swine fever virus and implications for viral assembly

A viral assembly factor promotes AAV2 capsid formation in the nucleolus

Viral assembly: a molecular modeling perspective

Interrogating viral capsid assembly with ion mobility-mass spectrometry

Coming together during viral assembly

Probing viral assembly

Synthesis, assembly, and processing of viral proteins

Thermodynamic basis for the genome to capsid charge relationship in viral encapsidation

Antibody Dynamics Simulation-Theory and Application

The matlab ode suite

SARS-CoV-2 D614G variant exhibits efficient replication ex vivo and transmission in 627 vivo

Spike mutation D614G alters SARS-CoV-2 fitness

Infection with the SARS-CoV-2 delta variant is associated with higher 630 infectious virus loads compared to the alpha variant in both unvaccinated and vaccinated individuals

SARS-CoV-2 Omicron variant shows less efficient replication and fusion activity when compared 632 with delta variant in TMPRSS2-expressed cells: Omicron variant replication kinetics

Novel insights into immune systems of bats

Weak protein− protein interactions are sufficient to drive assembly of hepatitis B virus capsids

Electrostatic repulsion, compensatory mutations, and long-range non-additive effects at the 639 dimerization interface of the HIV capsid protein

Physical regulation of the self-assembly of tobacco mosaic virus coat protein

Deciphering the kinetic mechanism of spontaneous self-assembly of icosahedral 646 capsids

Model-based analysis of assembly kinetics for virus capsids or other spherical polymers

Understanding the concentration dependence of viral capsid assembly kinetics-the origin of the lag 650 time and identifying the critical nucleus size

Physical principles in the construction of regular viruses[C]//Cold Spring Harbor symposia on 652 quantitative biology

Imaging and quantitation of a succession of transient intermediates reveal the 654 reversible self-assembly pathway of a simple icosahedral virus capsid

Genome-scale metabolic model of infection with SARS-CoV-2 mutants confirms guanylate 657 kinase as robust potential antiviral target

Molecular architecture of the SARS-CoV-2 virus

To bring data transparency in the era of COVID-19

Statistical Analysis Supports UTR (untranslated region) Deletion Theory in SARS-CoV-2

Emerging SARS-CoV-2 variants of concern (VOCs): An impending global crisis

Omicron sub-lineage BA. 2 may have "substantial growth advantage

Transmission of SARS-CoV-2 Omicron VOC subvariants BA. 1 and BA. 2: 666 Evidence from Danish Households