Impairment of T cells' antiviral and anti-inflammation immunities dominates the death from COVID-19

Clarifying key factors dominating the immune heterogeneity from non-survivors to survivors is crucial for therapeutics and vaccine developments against COVID-19. The main difficulty is to quantitatively analyze the multi-level clinical data of viral dynamics, immune response, and tissue damages. Here, we adopt top-down modeling to quantify key functional aspects and their dynamical interplays in the virus-immune system battle, yielding an accurate description of real-time clinical data involving hundreds of patients for the first time. The quantification of antiviral responses demonstrates T cells' dominant role in the virus clearance relative to antibodies, especially for mild patients(96.5%). Moreover, the anti-inflammatory responses, namely cytokine inhibition rate and tissue repair rate also have positive correlations with T cell number, and are significantly suppressed in non-survivors. Simulations show that impaired immune functions of T cells leads to greater inflammation (thus dominates the death), explaining the monotonous increase of COVID-19 mortality with age and higher mortality for males. We conclude that T cells play the role of crucial immunity that saves the death from COVID-19, which points out a new direction to advance current prevention and treatment by incorporating the vaccines, drugs and health care activities that aim to improve T cells' number and functions.


Introduction
The ongoing COVID-19 pandemic has resulted in over two million deaths worldwide. Therefore, identifying key factors that determine the immune heterogeneity from nonsurvivors to survivors is crucial for the current fight against the pandemic. Past clinical study has found a series of host factors associated with severe disease or higher mortality via correlation analysis: individual characteristics including older age, male sex, and comorbidities 1,2 ; profound lymphopenia, with T cells most significantly affected [3][4][5] ; the elevated level of inflammation markers, like LDH (lactate dehydrogenase) and Ddimer 2,6 ; excessive release of pro-inflammatory signaling molecules, like IFN − γ, IL-6, etc., known as the cytokine storm which is thought likely to be a major cause of multiorgan failure 4,7 . For immune responses, both SARS-CoV-2 specific T cells and antibodies are observed in COVID-19 patients 6,8 . However, the quantitative role of these factors in antiviral and anti-inflammatory immune responses is unknown, resulting in several unsolved questions about the cause and saving of the death: 1. What is the relative importance of T cell and antibody response for antiviral immunity at different stages; 2. What are the main driver and suppressors for the cytokine storm and multiorgan failure? And most importantly, 3.Are there new directions to overcome the heterogeneity of patients, decay of antibody function, and gene mutation of SARS-CoV-2 in efficient therapeutics and vaccine developments?
Beyond correlative analyses, quantitative modeling is a powerful tool to simulate the measured dynamical immune response to reveal the relative importance of different components 9 . However, most recent studies focus on the simple viral dynamics and its interactions with immune responses 10,11 and antiviral drugs [12][13][14][15][16][17][18][19][20][21] , without considerations of organ damages and disease progression. On the other hand, some multiscale simulations [22][23][24] incorporate existing knowledge about the viral dynamics, immune responses, and comorbidities to simulate the clinical outcomes. However, these approaches involve hundreds of model parameters, which have considerable value uncertainties that limit the reliability of predictions and systematic comparisons with clinical data. Therefore, previous studies are either oversimplified or overcomplicated to clarify the key factors dominating death.
In this work, we adopt a top-down modeling to quantify crucial aspects in the virusimmune system battle to overcome this difficulty. Here, the battle is classified into three kinds of functional behaviors, namely, the pathogenic functions (e.g., virus and inflammation), the protective functions (e.g., innate and adaptive immunity), and the organ damages. Integrating with the existing clinical and immunological knowledge for COVID-19 patients, we establish a dynamical motif for a small set of crucial functional variables and their interplays. This antiviral-inflammation model is used to simulate the systematic progression of COVID-19 patient with 19 measurable parameters. These simulations are validated with real-time clinical data involving hundreds of patients and then used to evaluate contributions of T cells and antibodies to antiviral immune responses. Subsequently, we quantify the difference of anti-inflammatory immune responses from non-survivors to survivors and clarify their correlations to T cells. Finally, T cells' dominant role in saving the death of COVID-19 and revelation to new therapeutics and vaccine development are discussed.

Causal network of the Antiviral-Inflammation Model
The difficulty of previous multiscale simulations 22-24 due to considerable parameter value uncertainties stems from the fact that, in the bottom-up strategy, the immune response to infectious disease is modeled as a complex network of numerous factors, resulting in the so-called 'curse of dimensionality' 25 . In contrast, a recent successful model of a classical complex system, namely, fluid turbulence, one of us has demonstrated that the global motions composed of numerous components typically display a symmetry-breaking which can be quantitatively modeled with finite functional variables, called order functions 26,27 . Here, we use a similar top-down strategy to quantify immune functional aspects in the virus-immune system battle. The model explicitly describes dynamics of five crucial functional quantities that determine COVID-19 progression: virus and interleukin 6 (IL-6) for pathogenic function, effector T cells and neutralizing antibodies for protective function, D-dimer (coagulation marker), and high-sensitivity cardiac troponin I (HSTC, heart injury marker) as examples for multi-organ damage. Other secondary factors modulate their interactions. Fig.1 shows their interplay following time order from virus dynamics, immune response to inflammation response. Self-replicating virus stimulates innate and adaptive immune cells which can produce antibodies. Effector T cells and neutralizing antibodies (NAbs) clear virus directly, either by killing infected cells or block the virus from entering into tissue cells; non-neutralizing antibodies (Non-NAbs) combines with innate immune cells (macrophages (Mɸ) and nature killers (NK), etc.) and induce immunoreaction to clear virus indirectly. The activated immune cells secrets cytokines, in which IL-6 has a central role for downstream destructive effect on organs, hence, it increases D-dimer and highsensitivity cardiac troponin I that characterizes multiorgan failure (MOF). On the other hand, suppressing immune hyperactivation and tissue repair by negative feedback . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint reduces activated immune cells and vessel and heart damage, thus reducing IL-6, Ddimer, and high-sensitivity cardiac troponin I.

Mathematical description of the Antiviral-Inflammation Model
Quantitatively, a set of ordinary differential equations, (1)-(5) is constructed to describe the antiviral immune response and pathogenesis of inflammation at a specific part of the body. The detailed meanings and units of all parameters are provided in SI (Supplementary Information), Table S1. The antiviral process includes continuous virial replication, dynamical activation of T cells to effector T cells by the virus, dynamical secretion of neutralizing antibodies by B cells, virus clearance by effector T cells and antibodies, and decay of effector T cells and antibodies. Eq. (1)-(3) describes the corresponding dynamical evolutions of concentrations of the virus, effector T cells, and neutralizing antibodies, with concise rate constants and time interval describing global effects of a series of microscopic processes, like target cell infection, antigen identification, antigen presentation, differentiation of T cells, B cell immunoglobulin class switching, etc. The antiviral effect of innate immunity is absorbed into an effective viral replication rate, (. Eq. (4) shows that the release of IL-6 is proportional to concentrations of non-neutralizing antibodies (assumed to be proportional to neutralizing antibodies) which combine with innate immune cells and trigger pro-inflammatory pathways. Immune suppressing reduces activated immune cells and, therefore, the IL-6 level. The release of IL-6 causes damage to organs, like thrombus formation and heart injury marked by D-dimer and high-sensitivity cardiac troponin I, whose levels are reduced when blood vessels and heart are repaired, shown in Eq. (5).
Simulations of Eq. (1)-(5) are compared to real-time data with 457 patients involved. Ten are individuals, and 447 patients form mild, severe, survivor, and non-survivor groups . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint whose median data are compared. For detailed data sources and integration of data from different sources, see Methods and SI. The time axis for simulation is the number of days after symptom onset, with the starting day being 0 or several days earlier. The numerical simulation is propagated using the delayed differential equations (DDE) of MATLAB. The least-square fit is achieved using the fmincon function of MATLAB with the implemented interior-point optimization algorithm. For a detailed description of fitting procedure and parameter uncertainty estimation, see Methods and SI.

Viral dynamics and contributions of T cells and antibodies to antiviral responses
Virus, effector T cell, and antibody dynamics are simulated and compared with real-time data from 10 individuals 28,29 and median values of mild and severe (critical) groups, survivors, and non-survivors [30][31][32] . The concentrations of virus (V), effector T cells (Te), and antibodies (A) is assumed proportional to viral load measurement from the respiratory tract, the change of T cells (T) from its baseline level (T0) measured in peripheral blood (T0 -T), and optical density or titer of Anti-RBD IgG/Anti-S1 IgG/Anti-NP IgG, respectively. CD3+ T cell (all T cells) data from cases of the same severity compensate for those who lack CD3+ T cell data. Original lymphocyte count is scaled by 0.589(Methods) for T cell count. As shown in Fig.2 and Supplementary Fig.1, simulations show agreement to data for all cases; for parameters, see SI.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint  Table S4. j: An overall statistic of the fraction of virus killed by T cells for all cases shows T cells' contribution is significantly higher in mild patients than severe patients. Solid markers are individual data and hollow markers are group data. Error bars represent standard errors. Clarifying deterministic factor controlling viral load peak benefits early antiviral treatment, vaccination, and epidemiological control 33 . By asymptotic analysis (Methods), the model predicts the peak is determined from virus inhibition by T cells: peak value is the ratio between square of virus replication rate and two times the multiplication of T cell activation and virus clearance rate. This gives 1.37 × 10 % copy/mL for the viral peak of patient P1 in Fig.2. Data agreement gives that for patients who survive, on average, 95% of viruses are cleared per day with 10 9 /L T cells from blood engaged, revealing strong efficiency of T cells' virus clearance. In conclusion, the consistency between simulation and data clarifies the antiviral dynamics for various severities in which adaptive response plays a significant role---first, effector T cells are activated, kill For patients of different severities, we compare the quantitative contributions of T cells and neutralizing antibodies for virus clearance as displayed in Fig.2j

Inflammation dynamics associated with death
To clarify the main driver for the cytokine storm and organ damage of critical illness, we compare simulations of Eq. (4)-(5) with real-time, median data of survivors and nonsurvivors (Fig.4). The concentration of non-neutralizing antibodies is assumed proportional to anti-RBD IgG optical density. For data source and parameter estimation, see Methods and SI. The agreement between simulations and data of IL-6, D-dimer, and high sensitive cardiac troponin I proves non-neutralizing antibodies stimulate the secretion of IL-6, and IL-6 stimulates the accumulation of D-dimer and high sensitive cardiac troponin I. IL-6 formation rates are assumed to be the same for both groups. The striking feature of non-survivors compared to survivors is the continuous production of IL-6 and organ damage, revealed by zero inhibition rates for all three markers (Fig.4b), while the difference of formation rates of organ damage markers is not remarkable.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint Figure 3. Comparison of predictions to clinical data of survivors and non-survivors reveals the mechanism for death from COVID-19. The data are the median of the group with error bars from reference 32 . The IL-6, D-dimer, and elevated sensitive cardiac troponin I inhibition rates are all zero for non-survivors, which gives continuous deterioration in contrast to non-zero inhibition rates of survivors who recover. For parameter estimation, see Methods and SI. The saturation value of Anti-RBD IgG for non-survivors is estimated by assuming its ratio to the maximum of survivors'(18 noncritical and three critical) data is close to the ratio of maximums of neutralizing antibodies between critical and non-critical patients 34 .
We conclude that the main trigger of cytokine storm is the secretion of elevated nonneutralizing antibodies which combine with innate immune cells and activate cytokine release pathways; the lack of negative feedback for immune suppression in non-survivors makes cytokine keeps accumulating. The cytokines act as the primary source that causes damage to organs, and the lack of tissue repair in non-survivors continues to worsen the illness. Our finding reveals the crucial aspect of death from COVID-19 is the lack of negative feedback for anti-inflammatory cytokine inhibition and tissue repair.

Initial T cells as background immunity that reduce mortality
Our model reveals that T cells' virus clearance, cytokine inhibition, and tissue repair are three essential protective functions in COVID-19 that determine disease severity. To seek what determines these protective functions, T cells' contribution of total virus clearance, cytokine inhibition rate, and tissue repair rates are plotted with initial T cell concentration before infection (equal to T cell baseline, T0) in Fig.4.a, is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint correlation (Fig. b is based on statistics of 137 non-survivors and 54 survivors). It shows the importance of sufficient initial T cells for comprehensive protection, which comes from an adequate number of effective CD8+ effector T cells that kill infected cells, sufficient regulatory T cells and other subsets that suppress the immune response and promote tissue repair 35 to reduce over inflammation in non-survivors. A great public concern is how an individual patient's background 'immune health' landscape (simplified as background immunity) shapes responses to SARS-CoV-2 infection 36 and controls the disease's severity. Because of the determining relation of initial T cell concentration for protective functions, T cells' static reserve before infection and dominance in population compared to other cells against the virus, we propose is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint concentration of initial T cells is a crucial characterization for the background immunity against SARS-CoV-2. To verify this hypothesis, we conduct disease progression of patients with different initial T cell concentrations. According to Fig.4 a,b, by assuming a linear decreasing of T cells' virus clearing rate, IL-6 and D-dimer inhibition rates with decreasing initial T cell concentration, Fig. 4c shows coagulation becomes more and more significant, which means lack of initial T cells exacerbates disease severity and increase mortality risk.
Therefore, we conclude that the T cells' impaired antiviral and anti-inflammation functions are the main immune origin of death from COVID-19: the extremely low level of initial T cells in non-survivors results in weak antiviral, cytokine inhibition and tissue repair abilities as well as low tissue repair function; then it calls the elevated antibodies for compensation; as a result, the concomitant large amount of non-neutralizing antibodies amplifies the cytokine storm, leading to continued damage. Following this casual chain, according to the decrease of lymphocytes (hence decrease of T cells assuming T cell count proportional to lymphocyte count) with older age (Fig.4d), we predict straightforwardly older patients must have higher mortality than younger patients. Also, male patients should have higher mortality than female patients for their lower level of CD4+ T cells 37 . Our prediction explains the continuous increase of COVID-19 mortality with age and higher mortality for males 37 , shown in Fig. 4d.

Discussion
In conclusion, we have quantified the adaptive-immune-response heterogeneity from non-survivors to survivors of COVID-19, using a dynamical motif with 19 measurable parameters beyond the overcomplication of the previous multiscale model 24 . For the first time, this model provides an accurate description of real-time clinical data involving hundreds of patients, which then reliably clarifies T cells' dominant roles in the antiviral and anti-inflammatory immune responses. Beyond the previous correlation analysis for T cell scarcity and disease severity 8,36 , this work reveals the causal relation between death from COVID-19 and impaired T cell immunity. This discovery explains the high mortality of older man and points out a new direction for therapeutics and vaccine development.
The currently tested drugs target various pathogenesis levels, from antiviral to antiinflammatory drugs and antithrombotic agents 24 , etc. However, there is no proven effective therapeutics for COVID-19. One crucial challenge is the lack of broad applicability of these drugs to heterogeneous patients with various comorbidities, disease severities, and complications 38 . Our study points out a new direction will be increasing T cell number and functions by both drugs and health care activities, which may benefit . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint virus clearance, cytokine inhibition, and tissue repair simultaneously. Firstly, recent studies reported the curing effect of drugs to COVID-19 patients by increasing T cell number, e.g., recombinant human granulocyte colony-stimulating factor 39 JinHuaQingGanKeLi 40 . Therefore, we encourage further studies and applications in this direction. Second, for the recovery of COVID-19 patients and healthy people's prevention, improving background immunity associated with T cells is more important and promising. Therefore, we strongly suggest studying the curing and immunity improvement effects of health care activities, such as mediation 41,42 , Tai Chi 43 and BaDuanJin 44 , for previous studies have found they help to increase CD3+ T cell and CD4+ T cell concentration and apply to a wide range of age, including older adults.
On the other hand, in the current development of the vaccine, neutralizing antibody immunity plays a crucial role. Unfortunately, the single-strand RNA structure makes SARS-CoV-2 easy to mutate, and several lineages have been discovered 45 . These mutations pose challenge for long-term effectiveness of antibody immunity, for it is on the molecular level targeting specific epitopes of the virus. By contrast, memory T cells show strong cross-reactivity and persistence 46 , and active T cells protect bodies in several aspects, including antivirus, suppress immune hyperactivation, and promote tissue repair. Therefore, stimulation of T cell response by the vaccine is worth more exploration, and we suggest advancing the current combination adjuvant strategy 47 that elicits potent CD8 and CD4 T cell responses.
For clinical application, to maximize the curing effect for severe patients, we suggest adopting multistage, synthetic protocols incorporating the above therapies. In this case, our model provides a strong tool to evaluate the effectiveness of treatments to identify individual optimal protocols. The reason is that all parameters can be determined from clinical data and quickly predict individual patients' trajectories, which may also advance the early prediction algorithm of current artificial intelligence softwares 48,49 .
To separate the critical elements from irrelevant details, we here have made some assumptions, which should be evaluated in further clinical studies, although they would not affect the basic conclusions of the study. For instance, the virus-clearance rate of innate immunity is thought to be a constant with negligible variation and is small compared to rates of adaptive response, which should be verified by further timedependent measurement of the course of innate response. The second assumption that needs more measurement to test is that, for survivors, the concentration of effector CD8+ T cells at the infected part is proportional to the reduction of CD3+ T cells in peripheral blood; whether it is strictly obeyed at the most time or not might be intriguing to further investigation 4 . The third questionable assumption is that the temporal profile of non-. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint neutralizing antibodies is proportional to neutralizing antibodies. Further measurement needs to test whether they are secreting at the same pace or not.

Methods § Asymptotic analysis of viral load peak
Around the viral load peak, with the decay of effector T cells and the production of antibodies neglected, the evolutions of viral load and effector T cells become: The analytical solutions are: All extracted data were made available to readers in our GitHub shared folder: https://github.com/luhaozhang/covid19 § Integration of data from different sources The data we used here contains 10 individuals and 4 groups. Individuals contain 6 mild and 4 severe (including critical) patients. The CD3+ T cell data of mild, severe and critical population from Hongzhou Lu 3 's cohort were used for individual patients and groups of the same severity who don't have T cell data. 4 patients were from Isabella Eckerle's 1 cohort and 6 patients from Kelvin To's 2 cohort, and the severity classification follows the assignments in previous publications. One patient from Kelvin To's cohort who have not yet been identified as being critically ill is classified as severe because the probability of being critically ill is low (1/6). For mild, severe (including critical) and non-survivor groups, Table S2 lists the sources of virus, T cells, and antibodies. T cells, IL-6, D-dimer and HSCT were from BinCao's cohort 32 for both survivors (137 patients) and non-survivors (54 patients). Virus and Anti-RBD IgG for survivors . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint were from the median of 20 survivors from Kelvin's cohort 50 . Lymphocyte data were multiplied by 0.589 to get estimated T cell concentration (0.589 is the ratio between medians of normal ranges of T cells 3

and lymphocytes 28 ) § Least square fit of virus, immune response and inflammation data
For parameters of simulations in Fig.2 and Fig.4, we adopt a best-fit approach to find the parameters which minimize the given objective function: the mean of residual sum of squares (RSM) between data points and the corresponding model simulations as used similarly in influenza model 51 . For virus-T cell-antibody dynamics, the objective function is:  (10) A is viral load data, B A is the value given by model, D is the total number of data points. The instructions are similar for T cells(T) and antibodies(A). For the objective function of inflammation response, mean of RSM (Eq. (11)) was used in linear scale for survivors and log scale (Eq. (12) (12) For simulation of viral load dynamics, when V<10 copy/mL, it is thought to be cleared thoroughly at one time without further evolution and is set to be 1 copy/mL. Table S3 lists the data points used for performing the best fit of each case. Viral data points that represent negative results of virus were abandoned unless they were important indicators of the ending of viral activity. The viral data points associated with the second viral load peak in mild and severe group were not used for any fits because the phenomenon is not observed as the common feature of individual patients. The decay of CD3+ T cells of 902 and 910 patients after 30 days were not used for fit at present. The decay was attributed as a deviation from averaged behavior on group level because of the large 95% CI of the data.
When performing best fit related with I, Sd, and Sh, parameters of virus, T cells and antibodies had been fixed. Estimation of uncertainty of parameters is carried out by performing a series of fits after the best fit for each case. For uncertainties of , , H , H , J , J , as shown in Fig.4, a set of fits was performed, each time by varying one data point of IL-6, D-dimer or HSCT to its upper or lower limit. The uncertainties for each fitted parameter were estimated as the minimum and maximum of each parameter among all fits. The best fits and uncertainties were summarized in Table S5. 95%CI is not used for uncertainty here because the number of data points is small. For methods and results of uncertainties of antiviral parameters, see SI.
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint fmincon function of MATLAB (MathWorks, version 2012 and higher) with the implemented interior point optimization algorithm was used to perform the fits. It requires constraints for the parameters to be optimized and an initial guess. An empirical fitting is performed for each patient and group to identify ranges of parameters' constraints for optimization, shown in Table  S6 and Table S8. Patients with the same type of data and same severity category are set to have similar ranges. Table S7 and Table S8 show the empirical values of parameters as initial guesses for best fits. Random initial guess is not suitable here because the fit is sensitive to the initial space, probably because of the limited number of data points with relatively large fluctuation, especially for viral load. For instructions about fixed parameters during the fitting, see SI. Parameters of best fits were used as initial guesses of fits when estimating parameters uncertainties. § Simulation of D-dimer dynamics with different initial T cell concentration Based on the positive correlation of T cell antiviral contribution, IL-6 and D-dimer inhibition rates, we assume a linear increase of T cell activation rate, and inflammation inhibition rates, , K , L with initial T cell concentration T0 as shown in the eq.(13)- (16). When T0 takes value of survivors' and non-survivors' median as in Fig.3(b), the parameters equal to their fitted values. Simulation in Fig.4(c) is performed using a series of given T0, the corresponding , , K , L and all other parameters of non-survivor group and their values remain unchanged. . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted April 29, 2021. ; https://doi.org/10.1101/2021.04.26.21256093 doi: medRxiv preprint