Mathematical model for assessing the level of cross-immunity between strains of influenza virus subtype H₃N₂

Full Text

Introduction

The World Health Organization (WHO) Global Influenza Surveillance and Response Network (GISRS) monitors and analyzes the evolution and epidemiology of influenza viruses with the main goal of selecting a vaccine strain and improving this process through research aimed at a better understanding of evolutionary variability, distribution factors in combination with the immunological landscape of the population and cross-immunity [1]. At the same time, WHO regularly updates recommendations for influenza vaccines in order to achieve their maximum compliance with the next circulating strains.

Nevertheless, over several seasons, the effectiveness of the influenza A vaccine, namely its H₃N₂ component, was determined to be low compared to other strains [2–5]. The reason for the low efficiency can be several factors. For example, characteristic adaptive changes during passaging of the recommended strain in chicken embryos during the production of vaccines [6]. Since the development, large-scale production and distribution of a vaccine takes many months, eventually the prevalence of circulating strains in some seasons will change significantly by the coming season [7]. But even if the recommended vaccine is well matched against circulating strains, its effectiveness could be adversely affected by the existing immune landscape [ 8]. Therefore, prediction of the evolutionary variability of the influenza virus is still of great interest for public health [9–12].

The most promising direction in this area is the construction of computer models that can be used to combine various modeling approaches, use multiple data sources with the ability to interpret the results for recommendations when choosing a vaccine strain. This requires close cooperation between scientists from different fields and directions, working at all levels of epidemiological surveillance and selection of vaccine strains, as well as model developers, epidemiologists and clinicians [10, 13, 14].

The team of the National Research Center of Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya in 2020 developed and successfully registered (certificate of registration No. 2020617965 dated July 15, 2020) the computer program Influenza IDE – an epidemiological model (EM) across continents with a simplified model of cross-immunity and a constantly updated database (of various types and subtypes of the influenza virus) Influenza DB. An important feature of EM is the possibility of forming the zero immune landscape of the population and then, after step-by-step simulation of the spread of the influenza virus, obtaining an immune landscape on the first day of the next season. The epidemiological model of the spread of the influenza virus among the world population over several seasons, developed using the agent approach, is presented in the form of implemented models: population behavior model, model of infectious process and infection model (based on the immune response in the body of an individual agent (person), taking into account the immune memory and the cross-immunity model). The computer program is designed to integrate various cross-immunity models.

A significant number of works on the study of cross-immunity is based on the discovery of a relationship between hemagglutination-inhibition antibody (HAI) titers and differences in the genetic sequences of viruses. To find this relationship, researchers use various mathematical methods, including regression analysis [15–18]. At the same time, various functions from HAI titers are used as a measure of cross-immunity. So, in the work of F.M. Burnet and D. Lush [19] introduced a function designated as an indicator of vaccine effectiveness R_ij = c_ij / c_ii, where c_ij is the concentration (dilution) value of serum to virus i in HAI assay with virus j, c_ii is the concentration (dilution) value of serum i in HAI assay with original virus, where the concentration (dilution) of serum in HAI assay is the reciprocal of the titer of HAI assay. Researchers I. Archetti and F.L. Horsfall [20] introduced the geometric mean of the above ratios (R_ijR_ji)^1/2 as a measure of antigenic variability in their work. In later works, which laid the foundation for the currently widely used antigenic cartography, A. Lapedes and R. Farber [21], D.J. Smith et al. [22] demonstrated that it is possible to construct a low-dimensional form space in which antibodies and antigens are considered as points, and the distance between them is denoted as antigenic distance. As a measure of distance D_ij, the logarithm of base 2 RTGA titers and the log₂(R_ij), 0,5 × log₂(R_ijR_ji) values were used [15–17].

As a measure of the difference between sequences, either a simple Hamming distance or a function depending on it is usually used. The most frequently used model in the considered works is the consolidation of amino acid positions into antigenic sites [23–25]. In such cases, the Hamming distance between antigenic sites is considered [16–18]. But there are also complex cases when physicochemical factors of differences in amino acid positions, for example, glycosylation, act as a measure of the difference between antigenic sites [18, 26]. As materials for building models, limited datasets were used, including only reference strains, which reduces objectivity. In this paper, we present a cross-immunity regression model based on the entire available array of WHO-published HAI data, which will improve the accuracy of the model, increase its objectivity and testability.

The use of an adequate epidemiological model with a verified cross-immunity model will improve the process of selecting the necessary vaccine strains for a more successful fight against the influenza virus.

The study purpose is to develop a mathematical model of cross-immunity based on the available array of seasonal serological testing data published by WHO (HAI assay).

Materials and methods

When forming the Influenza DB data array, information was used from WHO’s published seasonal data on the results of HAI testing and data (sequences + accompanying information) from the GISAID (Global Initiative on Sharing All Influenza Data) platform. The computer program Influenza IDE [27] also provides an additional module Influenza DP, which is a universal database processor (GISAID, NCBI, etc.).

To study patterns in the cross-immunity model, multiple regression analysis (linear regression) was used using the least squares method (LSM) and the non-negative LSM method [28] to estimate the regression parameters, represented by formula (1):

Modified_titre = c0 + c1 × As1 + c2 × As2 + c3 × As3 + c4 × As4 + c5 × As5 + c6 × As6, (1)

where the modified HAI titers are represented as the function value, as arguments As1, As2, As3, As3, As4, As5, As6 are the values of the degree of homology of antigenic sites (As) A, B, C, D, E, F, respectively, and c0, c1, c2, c3, c4, c5, c6 – parameters (coefficients) of the model.

For regression arguments, we relied on our own method for determining antigenic sites. The term “antigenic site” was introduced by Gerhard and Webster in 1978 to describe specific monoclonal antibodies. Antibodies that competed with each other were considered to bind the same antigenic site. Each antigenic site may contain one or more epitopes – different sets of amino acids on the antigen that come into contact with the amino acids of the antibodies. Competition between antibodies that bind the same antigenic site suggests that the epitopes at a given site physically overlap, but may be different, and one antibody molecule shields the entire antigenic site.

Also, to account for differences between strains by antigenic sites, we introduced the concept of homology degree and defined it for two randomly selected strains as a decreasing function of the Hamming distance, taking values from 1 to 0. The justification and selection of the functions used are described in detail in the section “Results”. The Stats package version 4.0.3 of the R programming language was used to train the regression model. The Pandas package version 1.4.2 of the Python programming language was used for data preprocessing. The stability of the predictive ability of the cross-immunity model was analyzed on retrospective data. As a measure of the adequacy of the model, we used the coefficient of determination R₂ [29], and as a measure of the accuracy of the forecast – the index of reproducibility of titles for one (± 1) and (or) two (±2) dilutions. We relied on numerous studies in which, when comparing the results of HAI assay for each sample within one laboratory or between several, the titers were considered equivalent if they differed by no more than one dilution (i.e. by 2 times) [30–32].

Building an Influenza DB Dataset to Study a Cross-Immunity Model

Preparation of experimental HAI data

WHO reports on influenza have been published every season since 2005 and are in the public domain [33]. The presented antibody titers in HAI reveal the antigenic properties of reference and test strains based on the cross-reactivity of ferret antibodies against reference strains that inhibit the agglutination of guinea pig erythrocytes by the tested influenza viruses. Due to the fact that WHO in 2008-2009 found a noticeable effect of neuraminidase and various sources of erythrocytes (turkey, human, guinea pig) on the results of HAI assay, in further studies we used only HAI assay tables from 2009 to 2022 performed on guinea pig erythrocytes with the addition of 20nM Oseltamivir, which is used to exclude the influence of neuraminidase. Despite the proposed international coding of viruses, researchers from all over the world encounter inaccurate annotations when depositing sequences in public databases [34]. Therefore, a separate task in our work was the procedure for processing and preparing published HAI assay tables for further use. Also, as a result of the preparation, each sequence was assigned a unique identifier.

Modified HAI titers act as a measure of cross-immunity in our model. To do this, the values of the published HAI titers are converted to the logarithm to the base 2. Next, the data is averaged according to the algorithm: subsets are combined in accordance with similar pairs of strain identifiers (reference and test); in each obtained subset, the arithmetic mean of the modified titers is calculated for the same pairs of strain identifiers.

It should be noted that the cultivation (passage) of viruses performed in chicken embryos may cause characteristic adaptive substitutions that change the receptor specificity of viruses and the nature of glycosylation [6, 35] and, as a result, affect cross-reactivity in HAI data. Therefore, as a result of the analysis of the primary data, we decided to single out several subsets in the initial HAI data array:

– with passage history in chick embryos (Egg-Egg 4226);

– with passage history in cell culture (Cell-Cell 28 621);

– with mixed passage history (Egg_s-Cell_t 16 463): reference viruses (to obtain control antisera) cultured in chick embryos in combination with test viruses grown in cell culture;

– with mixed passage history (Cell_s-Egg_t 5032): reference viruses (to obtain control antisera) grown in cell culture in combination with test viruses cultured in chick embryos.

Our solution was based on numerous comparative studies of antigenic mutations of the H3N2 influenza virus during cultivation in chicken embryos and cell passages [36–38]. Considering that the amount of data with a passage history of Cell-Cell and Egg_s-Cell_t is several times larger than the dimension of subsets with passage histories of Egg-Egg and Cell_s-Egg_t, we performed regression analysis on the data of the first two subsets.

Preparing data from the GISAID platform

The GISAID international platform was launched in 2008 and has since offered a reliable mechanism for the exchange of all genetic and influenza-related data for researchers, scientists and healthcare professionals [39]. As in the case of HAI data, when depositing sequences (when manually entering data, etc.), the format of the downloaded data is often distorted and misclassified. This, in turn, makes it difficult to analyze and process the available information to improve approaches (methods) to the selection of vaccine strains. Therefore, for further research using the computer program Influenza IDE and the universal processor Influenza DP built into it, we carried out a thorough cleaning and reconciliation of the available data. After the completion of the data verification process, the amino acid sequences were aligned to the corresponding reference sequence: H₃N₂: A/Aichi/2/1968, 566 AA, including the signal peptide, using the original fast lexical algorithm (maximum amino acid match between the reference sequence and the sample with minimization of the number of positions deletions and insertions) or the Smith–Waterman method [40]. Samples with differences of more than 20% from the reference were discarded. Brief statistical characteristics of the data are presented in Table 1.

 

Table 1. Data statistical characteristic

Таблица 1. Статистическая характеристика данных

Parameter Показатель		Value Значение
		Cell-Cell	Egg_s-Cell_t
Total number of observations Общее количество наблюдений		28 621	16 463
Total number of pairs of strains Общее количество пар штаммов		21 394	12 444
Number of observations in the reporting period Количество наблюдений в отчетном периоде	2009	21	35
	2010	46	94
	2011	132	461
	2012	2142	1761
	2013	1738	2447
	2014	2891	2201
	2015	1907	1520
	2016	1209	1003
	2017	1948	1279
	2018	659	164
	2019	2236	950
	2020	2496	1064
	2021	902	281
	2022	10294	3203
Titer value, Ме [Q1; Q3] Значение титра		160 [40; 320]	160 [40; 320]
Distribution by titer values Распределение по значениям титра	20	3404 (11,9%)	1858 (11,3%)
	40	4154 (14,5%)	2514 (15,3%)
	80	4232 (14,8%)	3852 (23,4%)
	160	6328 (22,1%)	3794 (23%)
	320	6170 (21,6%)	2446 (14,9%)
	640	3491 (12,2%)	1410 (8,6%)
	1280	772 (2,7%)	511 (3,1%)
	2560	60 (0,2%)	75 (0,5%)
	5120	10 (< 0,1%)	3 (< 0,1%)
Hamming distance across antigenic sites Дистанция Хэмминга по антигенным сайтам	A	2 [0; 4]	2 [0; 3]
	B	2 [1; 4]	3 [2; 4]
	C	1 [0; 2]	1 [0; 2]
	D	1 [0; 1]	1 [1; 2]
	E	1 [0; 2]	1 [0; 1]
	F	1 [0; 2]	1 [0; 2]

 

Results

Determination of antigenic sites of H₃N₂ influenza virus hemagglutinin

The combination of amino acid sequences into antigenic sites was carried out according to our own method with the inclusion of an additional antigenic site F. To determine which antigenic site of the influenza virus a particular position belongs to, we primarily focused on experimental studies [24]. Subsequently, we expanded the number of positions in antigenic sites by examining the variability of each hemagglutinin position starting from 1968. In the array of sequences obtained at the first stage, the occurrence of each amino acid in each position was counted. In parallel, an analysis was carried out for the ability of one or another amino acid of hemagglutinin to interact with antibodies, determining the immersion or, conversely, the exposure of amino acids on the surface of hemagglutinin. During the work, the GETAREA software [41] was used, the input of which was the tertiary structures of hemagglutinin molecules from the PDB (Protein Data Bank) database [42]. As a result, based on the variability of positions, their exposure in the tertiary structure of the protein on the surface, and also taking into account their maximum proximity to existing antigenic sites, variable exposed positions were determined, which were assigned to existing antigenic sites. Moreover, in addition to the well-known antigenic sites A, B, C, D, and E, we identified another antigenic site, respectively, named F. This antigenic site consists of many similar epitopes in the stem part of the hemagglutinin molecule. The database of tertiary structures contains more than 20 hemagglutinin complexes with monoclonal antibodies to this site. This site is described in detail in a study conducted by D.N. Shcherbinin et al. [43].

Thus, 6 antigenic sites were identified for hemagglutinin subtype H3. Below is a list of amino acid positions that make up these antigenic sites, a total of 109 positions (Table 2). Numbering is given by mature hemagglutinin H₃ [44]. It should be noted that the antigenic site F, unlike sites A–E, is located mainly in the HA₂ subunit of hemagglutinin.

 

Table 2. Amino acid positions at antigenic sites

Таблица 2. Аминокислотные позиции в антигенных сайтах

Site Сайт	Amino acid position numbers Номера аминокислотных позиций
A	121, 122, 124, 126, 131, 133, 135, 137, 138, 140, 142, 143, 144, 145, 146
B	128, 155, 156, 157, 158, 159, 160, 163, 186, 188, 189, 190, 192, 193, 194, 196, 197, 198, 199
C	45, 47, 48, 49, 50, 53, 54, 271, 273, 275, 276, 278, 280, 304, 311, 312
D	171, 172, 173, 201, 202, 207, 208, 212, 213, 214, 217, 219, 220, 221, 222, 223, 225, 226, 227, 229, 230, 242, 244, 246, 248
E	57, 58, 62, 63, 75, 78, 82, 83, 88, 92, 94, 96, 260, 261, 262
F	25, 31, 33, 326, 347, 361, 375, 384, 386, 450, 452, 453, 476, 479, 484, 487, 489, 501, 505

 

The information obtained was uploaded to the Influenza DB database and served as a template for combining the existing array with aligned amino acid sequences into antigenic sites. Further, the array with antigenic sequences was expanded by adding information about the content of the number of amino acid substitutions between any two strains in each of the 6 antigenic sites. As a result of the above actions, we obtained for each of the 6 antigenic sites for randomly selected or separate (each) pair of strains both the value of HAI titers and the number of amino acid substitutions (Hamming distance). At the same time, substitutions in antigenic sites, regardless of the type and specific position, are accepted to be equivalent. In our model, sequences that do not have substitutions in antigenic sites, but have changes in other positions, were considered antigenically identical.

Development of a mathematical model of cross-immunity of the influenza virus

Analysis and selection of a function to assess the degree of homology

The subsets with the largest number of observations were selected for regression analysis: Cell-Cell 28,621 and Egg_s-Cell_t 16,463. Also, instead of a simple Hamming distance, we decided to use a function from it, thereby introducing the concept of the degree of homology between two arbitrary antigenic sequences, indicating how close the antigenic sites of two different strains are to each other in terms of antigenic properties, and taking the value of this function. It was necessary to choose a function that could meet the requirements and would not contradict biological processes: in a situation of complete (absolute) homology of two strains (in the absence of substitutions in antigenic sites), the function takes the value 1, and in the absence of homology (changes were noted in each amino acid position included in the antigenic site) the value of the function must be equal to 0 (or close to 0). Obviously, the function should decrease with increasing Hamming distance.

In addition, as shown in [45], the function that depends on the Hamming distance and describes the cross-reactivity of antigens is concave. Based on this, we introduced an additional restriction for the function we are considering, namely, that it is concave or at least not convex, i.e. each subsequent substitution contributes no more to the decrease in homology than the previous one.

According to the results of studies by colleagues [16, 45], which showed that when the number of amino acid substitutions in the antigenic site is equal to or greater than 7, cross-reactivity of recognizing antibodies between strains is practically not observed, we introduced an additional requirement for the desired function: for 7 substitutions, the value function should take such a small value that it can be neglected. We also used the function 1 – x / 8 as the top-limiting function under consideration (for the values x = {1; 7}), and we took the value 1 – 7 / 8 = 0.125 as the significance threshold for 7 substitutions.

Taking into account the specified requirements, the following functions were considered (Figure):

• 1 – х / 8;
• exp(–a × x), where a = 1/3; 1/2; 1;
• 1 / (a × x + 1), where a = 1; 2; 3.

 

Figure. Functions for evaluating the degree of homology.

 

The values of HAI titers indicated in the original WHO tables as < (<40) in our model were replaced by a value of 20. We take a fixed value of the smallest value of the modified titer as the c0 coefficient, i.e. log2(20) ≈ 4.322.

The left part of Table 3 presents the results for a subset with a passage history on cell culture (Cell-Cell); in the right part of the table – with a mixed passage history (Egg_s-Cell_t).

 

Table 3. R² values by passage histories for the considered functions

Таблица 3. Значения R² по пассажным историям для рассматриваемых функций

Function Функция	Cell-Cell (21 580 pairs of strains 2009–2022) (21 580 пар штаммов 2009–2022)						Egg_s-Cell_t (12 444 pairs of strains 2009–2022) (12 444 пары штаммов 2009–2022)
	1_20%	2_20%	3_20%	4_20%	5_20%	ср._20%	1_20%	2_20%	3_20%	4_20%	5_20%	ср._20%
1 – x / 8	0,457	0,483	0,471	0,486	0,490	0,477	0,273	0,252	0,277	0,282	0,255	0,268
exp(–x / 3)	0,494	0,519	0,504	0,520	0,530	0,513	0,326	0,293	0,317	0,336	0,295	0,313
exp(–x / 2)	0,476	0,498	0,478	0,494	0,509	0,491	0,331	0,294	0,313	0,328	0,293	0,312
exp(–x)	0,366	0,380	0,346	0,363	0,393	0,370	0,253	0,205	0,200	0,242	0,207	0,222
1 / (x + 1)	0,441	0,461	0,438	0,453	0,471	0,453	0,312	0,272	0,283	0,315	0,275	0,292
1 / (2x + 1)	0,376	0,391	0,358	0,375	0,399	0,380	0,268	0,220	0,216	0,258	0,224	0,237
1 / (3x + 1)	0,328	0,340	0,301	0,320	0,348	0,327	0,213	0,162	0,149	0,196	0,168	0,178

 

Consideration of the selected functions was carried out by cross-validation: the entire amount of data was randomly divided 5 times in the ratio of 80% to 20%, where the first part was taken as a training sample, and the second – for a test one. The results were evaluated by the highest average value of the adjusted R² obtained for the test samples. We also ranked the results for all functions. When comparing the results between subsets of Cell-Cell and Egg_s-Cell_t, almost complete correspondence of ranks (Spearman’s rank correlation coefficient ρ = 0.96) of functions was observed in relation to subsets. Meanwhile, when comparing the values of the coefficients of determination for two subsets, the corresponding highest values were noted for the Cell-Cell subset (by 1.6 times), which confirms the expediency of separating the available data by passages, and also confirms that passages on chicken embryos contain characteristic adaptive substitutions, changing the receptor specificity of viruses and the nature of glycosylation [6, 35] and, as a result, affect cross-reactivity in HAI assay. Therefore, we carried out further calculations on the Cell-Cell subset. As a function describing the degree of homology, exp (–x / 3) was chosen, which showed the best result for the Cell-Cell subset and the best result for the Egg_s-Cell_t subset.

The results of regression analysis with the Cell-Cell subset (21,580 observations after averaging over pairs of strains) using the selected function with the requirement that the coefficients be non-negative are presented in Table. 4. In the first column of the table, in addition to the corrected R² already mentioned above, the designations of the corresponding antigenic sites are presented, then in the Estimate column the values of the desired C1–6 regression coefficients are given, in the third column, Standard error, the standard deviations of these coefficients are listed. The line Standard deviation represents the standard deviation of the dependent variable (modified titers). The contribution of each antigenic site can be judged by the value of its coefficient. The greater the value of the coefficient at a certain antigenic site, the greater its contribution to cross-immunity. After determining the coefficients of the regression model, the statistical hypothesis is tested about the equality to zero of the true values of the coefficients according to Student’s criterion with a significance level of 0.05. The results of the test are presented in the fifth column, Pr(> |t|) p-value, the probability that the t-statistic is greater than the t-value modulus (ratio of coefficient values to their standard deviation). If this value is less than the confidence level of 0.05, the hypothesis is rejected and the parameter is considered significant. In our case, all antigenic sites, except for C and F, were recognized as significant. The coefficients at the antigenic site C and F took the value 0, and this is equivalent to insignificance.

 

Table 4. Regression results for the exp(–x / 3) function

Таблица 4. Результаты регрессии для функции exp(–x / 3)

Аntigenic sites Антигенные сайты	Cell-Cell_2009–2022 (21 580 pairs of strains) (21 580 пар штаммов)
	estimate значение	standard error cтандартная ошибка	P-r(> |t|) p-value
А	1,233	0,034	<0,001
B	2,506	0,034	<0,001
C	0,000
D	0,083	0,027	0,002
E	0,622	0,039	<0,001
F	0,000
CR2	0,514
Standard deviation Стандартное отклонение	1,161

 

Analysis of the stability of the predictive ability of the cross-immunity model based on retrospective data

For a more detailed understanding of the influence of the time factor and the number of observations involved in tuning the model on the accuracy of the forecast, the next step is to divide the subset of Cell-Cell_2009–2022 into different time periods (1-, 2-, 3-,4-, 5-year) adjustable periods with a forecast for each subsequent year. The full results of the regression analysis are presented in Appendix.

Based on presented in table 5.1 and 5.2 of the results, we can state that the adjusted R2, which takes a value not higher than 0.3 in adjustable periods, begins to increase significantly from 2019, changing from 0.32 to 0.64 (the more years in the adjustable interval, the lower the value). 2021 is out of this trend due to low incidence and little data due to the COVID-19 pandemic. The dynamics of the value of the adjusted R² is also directly reflected in the forecast periods, taking values no higher than 0.4 (including negative ones) and further increasing from 2020 (from 0.47 to 0.6).

 

Table 5.1. R²

Таблица 5.1. R²

Year Год	1 year год	2 year года	3 year года	4 year года	5 year лет
2014	0,135	0,099	0,099	0,101	0,101
2015	–0,223	–0,204	–0,354	–0,355	–0,348
2016	0,043	–0,225	–0,278	–0,540	–0,546
2017	0,005	0,198	0,158	0,083	0,071
2018	0,092	0,001	0,100	0,071	0,052
2019	0,124	0,244	0,405	0,318	0,277
2020	0,596	0,592	0,577	0,561	0,543
2021	0,044	0,095	0,125	0,135	0,128
2022	0,600	0,472	0,497	0,540	0,561

 

Table 5.2. ±1 dilution

Таблица 5.2. ±1 разведение

Year Год	1 year год	2 year года	3 year года	4 year года	5 year лет
2014	67,84	67,94	67,94	67,24	67,24
2015	46,46	46,69	45,67	45,59	45,75
2016	40,43	38,41	36,50	34,48	34,36
2017	56,40	59,40	63,40	55,83	59,54
2018	65,27	59,62	62,97	60,88	61,09
2019	45,44	50,94	59,41	57,20	57,57
2020	57,70	57,03	58,33	57,32	56,79
2021	46,03	48,12	50,63	53,35	51,88
2022	55,27	50,18	54,16	55,66	60,49

 

We also evaluated the predictive ability of the results presented in Table. 5.1 and 5.2, according to the reproducibility of titers for one (±1) and two (±2) dilutions. When comparing the calculated titles with those published in the forecast periods, with the exception of 2016 and 2021 (a small number of observations – 841 and 478, respectively) and 2015 (when the recommended vaccine strains failed to provide effective protection), the reproducibility of titers per (±1) dilution was demonstrated in more than 50% of cases, and at times reached 67% of cases. As for the reproducibility of titers at two (±2) dilutions, approximately 90% of cases and above are consistently observed over all time periods.

The results of regression analysis over adjustable (1-5year time intervals) demonstrated the significance of the contribution of each antigenic site to the immune response. The coefficients at antigenic sites A and B have high values in all forecast periods. This can be seen especially clearly on 3–5-year customizable periods. Antigenic site C actually takes on null values only with rare exceptions. Almost at all time intervals, a stable contribution, commensurate with the contribution of A and B, also demonstrates the antigenic site D, but only until 2018, followed by a sharp zeroing. Moderate compared to antigenic sites A and B, but at the same time stable significance is demonstrated by the antigenic site E with a noticeable increase by 2018–2021. The contribution of the antigenic site F varies diametrically depending on time intervals, from stable significant to zero. It should be noted that, in general, the results of the regression do not reveal a significant dependence on the number of years included in the adjusted period. At the same time, forecast results are significantly worse in those seasons when there is little data in the forecast period. In addition, our correlation analysis showed that there is no correlation between the values of the coefficients at antigenic sites and the number of substitutions in positions. Taking into account the fact that the number of observations only for 2022 is more than 1/3 of the total volume of observations for 2009–2022, it was decided to conduct separate calculations for observations in 2022. Figure 6 shows the results of calculating the predictive ability of the developed cross-immunity model based on the published data of the RTGA 2022 (spring and autumn seasons).

 

Table 6. Regression results for 2022 data

Таблица 6. Результаты регрессии по данным 2022 г.

Customizable period Настраиваемые периоды	2022 г.	February 2022 report data Данные февральского отчета 2022 г.
Number of observations Количество наблюдений	8183	1994
Adjusted R2 Скорректированный R2	0,728	0,625
Standard deviation Стандартное отклонение	0,874	0,966
Forecast periods Прогнозные периоды		September 2022 report data Данные сентябрьского отчета 2022 г.
Number of observations Количество наблюдений		6675
Adjusted R2 Cкорректированный R2		0,734
Standard deviation Стандартное отклонение		0,933
±1 dilution / разведение		69,33%
±2 dilution / разведения		98,73%

 

Discussion

The analysis of the scientific literature demonstrated that in most studies the assessment of cross-immunity was carried out on a limited data set, including only reference strains, which reduces the objectivity of studies and makes it difficult to verify the forecast. Carrying out calculations based on the entire available array of HAI data published by WHO made it possible to bypass these limitations. With the help of the built-in additional module Influenza DP, which is a universal database processor (GISAID, NCBI, etc.), we have developed an algorithm that allows us to generate databases in real time according to the tasks for further analysis and calculations. To study the patterns in the cross–immunity model, subsets with the largest number of observations – Cell-Cell and Egg_s-Cell_t – were considered. We also introduced the concept of the degree of homology between two arbitrary strains, which takes the value of a function depending on the Hamming distance.

When considering functions for choosing a variant describing the degree of homology, a significant dependence of the values of the coefficient of determination on the type of function was noted. While all functions were subjected to the same necessary conditions corresponding to biological processes, the resulting adjusted R² values fluctuated quite significantly and took values from 0.327 to 0.513 for the Cell-Cell subset and 0.178 to 0.313 for the Egg_s-Cell_t subset, respectively. Such a significant difference between the upper and lower limits of the indicator indicates that the use of a simple Hamming distance [16, 18] as a measure of homology can impose significant restrictions on the models used. The ranking of the results for all functions was also carried out. When comparing the results with the Cell-Cell and Egg_s-Cell_t subsets, almost complete correspondence of the function ranks by the value of the adjusted R² in relation to the subsets was observed. The results obtained may indicate a certain correspondence to the biological mechanisms of our approach to the choice of function.

Comparison of the adjusted R² values on the Cell-Cell and Egg_s-Cell_t subsets (0.523 and 0.313, respectively) demonstrates the validity of our decision to split the original dataset according to passage history. The difference in the adjusted R² also suggests that viruses cultured in chick embryos introduce a large error in the results of HAI assay.

In our calculations, we relied on the results of our own studies of combining amino acid positions into antigenic sites, described an additional F site, thereby expanding the number of antigenic sites to 6. The results of regression analysis on all Cell-Cell subsets showed that antigenic sites A and B made the greatest contribution to the immune response. They are moderately overtaken in stable values by the antigenic site E with a noticeable increase in 2018–2021. Antigenic site C practically does not contribute to the immune response in all considered intervals. Antigenic site D, which largely confidently demonstrates its significance at the beginning of the forecast periods, begins to lose ground in 2018–2021. The antigenic site F with its modestly unstable contribution to cross immunity turned out to be quite sensitive to time intervals.

In addition, we analyzed the stability of the predictive ability of our model for different time intervals (from 2009 to 2022). After dividing the Cell-Cell subset into 1 to 5 year time intervals, it was noted that the value of the adjusted R² in adjustable periods increases starting from 2019 reaching around 0.5 and above, except for 2021, which has very few observations (478) due to the COVID-19 pandemic. A similar picture is observed in the forecast periods. The pronounced high value of the adjusted R² for later data may be due to an increase in ongoing research (HAI assay) on cell culture and an improvement in their quality.

When comparing calculated values with published titers, reproducibility results per (±1) dilution are fairly stable across forecast years, with values greater than 50% (averaging around 60%), except for 2015, 2016, and 2021. At the same time, we see the main reason for the poorer results in a small amount of data for the forecast periods.

As shown in the results, the unprecedented number of observations for 2022 provided us with the opportunity to tune the model to the data at the beginning of 2022 (adjusted R² = 0.728), and check the results of the forecast for the autumn season. The prediction results (adjusted R² = 0.734) demonstrated titer reproducibility per (±1) dilution in 69.33% of cases and turned out to be significantly better than all previously considered periods, which can already be considered sufficient for practical application. The explanation for this result can be both the quality of the data and the forecast period limited by the season. In the work of T. Bedford et al. [15], who performed a regression analysis on the entire H₃N₂ HAI data set, obtained R²= 0.372, which is lower than our adjusted R² = 0.523 on the entire Cell-Cell subset. The results of colleagues could be affected by the lack of separation of data by passage history and time intervals. Also, the results may have been influenced by the use of a simple Hamming distance. It is important to recognize that, despite all the efforts to standardize the conduct of RTGA [32] and regardless of the chosen approach for studying the dependence of HAI titers on amino acid substitutions, the initially inherent high error of the HAI method (17%) [46] remains a factor that significantly affects the results.

Conclusions

Based on an epidemiological model and with a well-tuned cross-immunity model, we have developed a computer program that provides the ability to predict the most common strains, taking into account the influence of the immune landscape and select a vaccine strain for the upcoming season.

The ability to carry out calculations based on the entire available array of WHO-published seasonal HAI data, the observed improvement in results over recent years (starting in 2019), as well as a large volume of observations and the results of the forecast for 2022 (adjusted R² = 0.734; reproducibility of titers per one (± 1) dilution in 69.33% of cases) allow to assert that the proposed method can serve as a good tool for future forecasts with further study to confirm its stability.

In turn, the advantage of using the Influenza IDE computer program with a constantly updated database of various types and subtypes of the influenza virus in the future will allow reproducing the results obtained on other variants and subtypes of the influenza virus and thereby testing the developed model and expanding the spectrum of recommended vaccine strains. The experience and skills gained in determining antigenic sites in the future will make it possible to combine amino acid positions into antigenic sites according to several scenarios and making calculations already taking them into account. Further development of the cross-immunity model can be achieved by unilateral normalization of the titers of the entire array of test strains against reference strains, which will reduce the error associated with the individual characteristics of the experimental animals used. Also, the improvement of the model can be achieved by taking into account processes such as glycosylation.

About the authors

Marina N. Asatryan

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Author for correspondence.
Email: masatryan@gamaleya.org
ORCID iD: 0000-0001-6273-8615

PhD (Med.), senior researcher epidemiological cybernetics group of the Epidemiology Department

Russian Federation, 123098, Moscow

Boris I. Timofeev

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0001-7425-0457

PhD (Phys.-Mat.), senior researcher D.I. Ivanovsky Institute of Virology Division

Russian Federation, 123098, Moscow

Ilya S. Shmyr

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-8514-5174

researcher epidemiological cybernetics group of the Epidemiology Department

Russian Federation, 123098, Moscow

Karlen R. Khachatryan

National Research University Higher School of Economics

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-1934-532X

master's student

Russian Federation, 123458, Moscow

Dmitrii N. Shcherbinin

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-8518-1669

PhD (Biol.), researcher, Department of Genetics and Molecular Biology of Bacteria

Russian Federation, 123098, Moscow

Tatiana A. Timofeeva

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-8991-8525

PhD (Biol.), head of laboratory D.I. Ivanovsky Institute of Virology Division

Russian Federation, 123098, Moscow

Elita R. Gerasimuk

State University “Dubna”

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-7364-163X

PhD (Med.), Assoc. Prof.

Russian Federation, 141982, Dubna

Vaagn G. Agasaryan

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0009-0009-3824-7061

researcher epidemiological cybernetics group of the Epidemiology Department

Russian Federation, 123098, Moscow

Ivan F. Ershov

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-3333-5347

researcher epidemiological cybernetics group of the Epidemiology Department

Russian Federation, 123098, Moscow

Tatyana I. Shashkova

Artificial Intelligence Research Institute

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-8754-8727

PhD (Biol.), senior researcher Bioinformatics group

Russian Federation, 121170, Moscow

Olga L. Kardymon

Artificial Intelligence Research Institute

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-4827-8891

head of Bioinformatics research group

Russian Federation, 121170, Moscow

Nikita V. Ivanisenko

Artificial Intelligence Research Institute

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-0333-8117

PhD (Biol.), senior researcher Bioinformatics group

Russian Federation, 121170, Moscow

Tatyana A. Semenenko

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0002-6686-9011

D. Sci. (Med.), Prof., Full Member of RANS, head Department of Epidemiology

Russian Federation, 123098, Moscow

Boris S. Naroditsky

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0001-5522-8238

D. Sci. (Biol.), professor, Deputy Director for research D.I. Ivanovsky Institute of Virology Division

Russian Federation, 123098, Moscow

Denis Yu. Logunov

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org

D. Sci. (Biol.), Full Member of RAS, Deputy Director for research

Russian Federation, 123098, Moscow

Aleksander L. Gintsburg

National Research Center for Epidemiology and Microbiology named after Honorary Academician N.F. Gamaleya

Email: masatryan@gamaleya.org
ORCID iD: 0000-0003-1769-5059

D. Sci. (Biol.), Prof., Full Member of RAS, Director

Russian Federation, 123098, Moscow

References

Supplementary files