Computational study of the furin cleavage domain of SARS-CoV-2: delta binds strongest of extant variants

We demonstrate that AlphaFold and AlphaFold Multimer, implemented within the ColabFold suite, can accurately predict the structures of the furin enzyme with known six residue inhibitory peptides. Noting the similarity of the peptide inhibitors to polybasic furin cleavage domain insertion region of the SARS-CoV-2, which begins at P681, we implement this approach to study the wild type furin cleavage domain for the virus and several mutants. We introduce mutations in silico for alpha, omicron, and delta variants, for several sequences which have been rarely observed, for sequences which have not yet been observed, for other coronaviruses (NL63, OC43, HUK1a, HUK1b, MERS, and 229E), and for the H5N1 flu. We show that interfacial hydrogen bonds between the furin cleavage domain and furin are a good measure of binding strength that correlate well with endpoint binding free energy estimates, and conclude that among all candidate viral sequences studied, delta is near the very top binding strength within statistical accuracy. However, the binding strength of several rare sequences match delta within statistical accuracy. We find that the furin S1 pocket is optimized for binding arginine as opposed to lysine. This residue, typically at sequence position five, contains the most hydrogen bonds to the furin, and hydrogen bond count for just this residue shows a strong positive correlation with the overall hydrogen bond count. We demonstrate that the root mean square backbone C-alpha fluctuation of the first residue in the furin cleavage domain has a strong negative correlation with the interfacial hydrogen bond count. We show by considering the variation with the number of basic residues that the maximum mean number of interfacial hydrogen bonds expected is 15.7 at 4 basic residues.