By Luciano Abadias, Gissell Estrada-Rodriguez, and Ernesto Estrada
Since its onset in December 2019, the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has caused an outbreak of pulmonary disease that soon morphed into the COVID-19 global pandemic. No specific drug currently exists to treat SARS-CoV-2, which makes the situation even more critical. Drug repurposing, a process through which scientists identify new uses for previously approved or investigational drugs, could offer a rapid means of drug discovery to fight SARS-CoV-2 and the associated COVID-19 disease [3, 6]. This method is advantageous because it begins with drugs that were previously found to be sufficiently safe in humans; most preclinical tests are complete and the necessary technology for production is available. Drug repurposing starts with hypothesis generation, wherein researchers identify the candidate molecule that targets a viral protein.
SARS-CoV-2 Mpro as a Target
The life cycle of SARS-CoV-2 commences when the virus enters a human cell through receptor-mediated endocytosis [3]. Once inside the cell’s interior, the virus releases its RNA and proteins in the cytoplasm. The viral messenger RNA is then translated into two viral replicase polyproteins, which are cleaved by two viral proteases—the papin-like protease and the main protease (Mpro)—to create 16 non-structural proteins (NSP). The role of Mpro in NSP production is fundamental for virus proliferation in the human cell. It therefore constitutes an important target against SARS-CoV-2.
X-ray crystallography has resolved the three-dimensional (3D) structures of free SARS-CoV-2 Mpro, as well as Mpro that is complexed with some inhibitors. Mpro consists of two identical chains, each with 306 amino acids (AAs). A comparison of the SARS-CoV-2 Mpro with that of SARS-CoV (the virus that produced the 2002-2003 outbreak) revealed several similarities. They share 96 percent of their residue sequences, and superposition of their 3D structures revealed an almost perfect fit. But are they truly similar, or does the SARS-CoV-2 Mpro hide some “magic tricks” that make it more efficient in its damaging work?
SARS-CoV-2 Mpro as a Network
One can represent a protein as a network in the following way. Let the \(\textrm{C}_{\alpha}\) atom of each residue be represented as a node on a graph, and let us fix a cutoff radius \(r_{C}\), which represents an upper limit of separation between two residues in physical contact. Two nodes \(i\) and \(j\) are thus connected if \(r_{ij}\leq r_{C}\). This graph is called a protein residue network [2]. Using this network representation, Estrada studied the 3D structures of the SARS-CoV and SARS-CoV-2 Mpro [4]. To compare both proteases, he utilized several network-theoretic descriptors: edge density, degree heterogeneity, clustering coefficients, modularity, average shortest path distance, and average betweenness centrality. His findings confirmed the previous observations that both proteases are very similar, with less than a two percent difference in the aforementioned descriptors. Other descriptors, such as communicability indices, experience up to a 20 percent difference between the two proteases. They are based on the following matrix function of the network’s adjacency matrix \(A\):
\[G:=\sum_{k=0}^{\infty}\frac{A^{k}}{k!}=\exp\left(A\right),\tag1\] where \(\left(A^{k}\right)_{ij}\) counts the number of walks of length \(k\) between two nodes \(i\) and \(j\).
However, these indices heavily penalize relatively long walks, due to the fast-growing nature of \(k!\) Therefore, the possibility of exploring the proteases’ capacities to transmit effects through longer distances was an incognita. So Estrada examined a modification of communicability, which employs a double rather than single factorial [5]:
\[Z:=\sum_{k=0}^{\infty}\frac{A^{k}}{k!!}=\frac{1}{2}\left[\sqrt{2\pi}\textrm{erf}\left(\frac{A}{\sqrt{2}}\right)+2I\right]\textrm{exp}\left(\frac{A^{2}}{2}\right).\tag2\]
The results from this new matrix function were astonishing. Namely, it indicated that the SARS-CoV-2 Mpro is 1900 percent more likely to transmit long-range effects than its SARS-CoV analogue. In addition, the AAs displayed increased sensitivity to perturbations around the new protease’s binding site and near its catalytic site.
Figure 1. Schematic illustration of the main findings reported in [4]. 1a. The SARS-CoV-2 main protease (Mpro). 1b. The protein residue network of Mpro. 1c. Contribution of amino acids (AAs) near the binding site of the matrix function \(Z\). Figure courtesy of Ernesto Estrada.
Estrada used these findings to analyze several structures of the SARS-CoV-2 protease, which was bounded to inhibitors. Remarkably, the strongest inhibitors of the SARS-CoV-2 Mpro are those that produce the least amount of change in its capacity to transmit perturbations across the protein, thus identifying a potential route for researchers’ design of potent inhibitors for the SARS-CoV-2 Mpro. Figure 1 illustrates these results.
Networks-Fractional Calculus Alliance
A question still remained about the interactions of Mpro inhibitors: What is the mechanism through which one residue transmits its perturbations, such that it is “felt” at long distances in the protein? Although the communicability function \(G=\exp(A)\) can be directly related to dynamical processes on the graph, the same is not true of the function \(Z\). Thus, the question is actually as follows: Does a dynamical model exist on the network whose solution is related to \(Z\)?
The answer to this query stems from fractional calculus, a field that allows researchers to capture the sub- and superdiffusive behavior of many real-life processes. We developed the following fractional-order susceptible-infected (FSI) model to study the propagation of perturbations produced by protein inhibitors [1]:
\[\left\{ \begin{array}{l}
{D_{t}^{\alpha}(-\log(1-x))(t)=\beta^{\alpha}Ax(t),}\\
{x(0)=x_{0}.}
\end{array}\right.\]
Figure 2. Values of the average changes in transmissibility of perturbations \(\varDelta E^{1/2}\), average path lengths between the most perturbed amino acids (AAs) \(L\), and number of residues in these paths that belong to the binding site of the main protease (Mpro) \(N_{BS}\) for the interaction between Mpro and three different inhibitors. \(IC_{50}\left(\mu M\right)\) is the inhibitory concentration of the inhibitor. The least powerful inhibitor is 6Y2G and the most powerful is 6M0K.
Here, \(x(t)\) is a vector quantity that depends on time \(t\), which represents the probabilities that the corresponding residues will become perturbed at \(t\). \(\beta\) is a positive constant that signifies the rate at which the perturbation is transmitted between AAs, and \(\alpha\) is a number that belongs to the interval \((0,1)\), which provides the fractional order of derivation. We then employed the fractional Caputo derivative \(D_{t}^{\alpha}\). For example, \(\alpha=1\) yields the classical susceptible-infected (SI) model, which one can rewrite as \(x'(t)=-(1-x(t))\beta Ax(t)\).
We found an explicit upper bound \(\hat{x}\) to the FSI model’s exact solution, in terms of Mittag-Leffler functions of \(A\). That is, \(\hat{x}=f(\hat{y})\) where \(f(\hat{y})=1-e^{-\hat{y}}\) and
\[\hat{y}\left(t\right)=E_{\alpha,1}\left((\beta t)^{\alpha}\hat{A}\right)g\left(x_{0}\right)+\sum_{k=0}^{\infty}\frac{(\beta t)^{\alpha\left(k+1\right)}\hat{A}^{k}Ab\left(x_{0}\right)}{\Gamma\left(\alpha\left(k+1\right)+1\right)},\]
with \(E_{\alpha,1}\) as the Mittag-Leffler function, \(\hat{A}=A\text{diag}\left(1-x_{0}\right)\), \(b(x_{0})=x_{0}+(1-x_{0})\log(1-x_{0})\), and \(g=f^{-1}\). The upper bound \(\hat{x}\) represents the worst-case scenario on the propagation of perturbations across a protein residue network because \(0\leq x_{i}(t)\leq\hat{x}_{i}(t)\). Moreover, \(\hat{x}\) converges asymptotically to the exact solution as time approaches infinity. When the perturbation is produced with equal probability at any residue of the protein, this upper bound is written as
\[\hat{y}\left(t\right)=\nu E_{\alpha,1}\left(\zeta A\right)\vec{1}-\phi\vec{1},\tag3\]
where \(\nu\) and \(\phi\) are constants for given fixed parameters and \(\eta=(\beta t){}^{\alpha}\text{diag}(1-x_{0})\).
Figure 3. Illustration of the shortest paths between the most perturbed pairs of amino acids (AAs) in the main protease (Mpro), bounded to two inhibitors of different potency. 6Y2G is the least powerful inhibitor and 6M0K is the most powerful inhibitor. Figure courtesy of Ernesto Estrada.
We studied the average differences \(\varDelta E^{\alpha}\) between all nondiagonal entries of \(E_{\alpha,1}\left(\zeta A\right)\) when the protein is bounded to an inhibitor and when it is free for different values of \(\alpha\). For each \(\alpha\), we obtained the average length \(L\) of the paths between the pairs of nodes in the top ten ranking according to \(\varDelta E^{\alpha}\). We then counted the number \(N_{BS}\) of nodes in these rankings that belonged to the SARS-CoV-2 Mpro’s binding site. When \(\alpha=1\) (a normal SI model), the results indicated no clear trend between the calculated parameters and the potency of the inhibitors, as determined by \(IC_{50}\left(\mu M\right)\). But when \(\alpha=0.5\), the most powerful inhibitor increases the transmissibility of perturbations through Mpro by 71 percent after its binding. The second most powerful inhibitor increases the transmissibility of perturbations by 13 percent. Interestingly, the weakest inhibitor actually lessens the transmissibility of perturbations across the protein (see Figure 2). In addition, the average length of the shortest paths and number of residues in the binding site follow the same trend as \(IC_{50}\left(\mu M\right)\). Figure 3 offers an illustration of these perturbing paths for the least and most powerful inhibitors.
The resulting Mittag-Leffler function of the network’s adjacency matrix from the FSI model is
\[E_{1/2,1}\left(A\right)=\left[\textrm{erf}\left(A\right)+I\right]\textrm{exp}\left(A^{2}\right)=\sum_{k=0}^{\infty}\frac{A^{k}}{\Gamma\left(k/2+1\right)},\tag4\]
and the connection with the function \(Z\) stems from the fact that
\[k!!=2^{\left(1+2k-\cos\left(k\pi\right)\right)/4}\pi^{\left(\cos\left(k\pi\right)-1\right)/4}\Gamma\left(k/2+1\right).\tag5\]
These efforts collectively illustrate the use of applied mathematics in helping to understand and combat COVID-19.
References
[1] Abadias, L., Estrada-Rodriguez, G., & Estrada, E. (2020). Fractional-order susceptible-infected model: Definition and applications to the study of COVID-19 main protease. Fract. Calc. Appl. Anal., 23, 635-655.
[2] Estrada, E. (2011). Protein residue networks. In The Structure of Complex Networks: Theory and Applications (pp. 277-301). New York, NY: Oxford University Press.
[3] Estrada, E. (2020). COVID-19 and SARS-CoV-2: Modeling the present, looking at the future. Phys. Rep. To be published.
[4] Estrada, E. (2020). Topological analysis of SARS CoV-2 main protease. Chaos, 30, 061102.
[5] Estrada, E., & Silver, G. (2017). Accounting for the role of long walks on networks via a new matrix function. J. Math. Anal. Appl., 449, 1581-600.
[6] Nguyen, D.D., & Wei, G.-W., (2020, May 1). Math and AI-based repositioning of existing drugs for COVID-19. SIAM News, 53(4), 5.
Luciano Abadias is an associate professor in the Department of Mathematics and the Institute of Mathematics and Applications (IUMA) at the University of Zaragoza in Spain. His research is devoted to operator theory and applications, as well as fractional calculus. Gissell Estrada-Rodriguez is a postdoctoral research fellow in mathematics at Sorbonne Universite’s Laboratoire Jacques-Louis Lions in Paris. She works on the derivation and analysis of nonlocal partial differential equations with applications to biology and interacting particle systems. Ernesto Estrada is a 2019 SIAM Fellow and a senior researcher at the University of Zaragoza’s IUMA and the Government of Aragon’s ARAID (Agencia Aragonesa para la Investigación y Desarrollo) Foundation in Spain. He has a wide range of interests in the mathematics of discrete complex systems.