key: cord-0822877-yji7tj0w authors: Soni, Bhavnita; Singh, Shailza title: COVID-19 co-infection mathematical model as guided through signaling structural framework date: 2021-03-26 journal: Comput Struct Biotechnol J DOI: 10.1016/j.csbj.2021.03.028 sha: 4d6262c75a3bb30c2e448c57eb7c37ccb5001621 doc_id: 822877 cord_uid: yji7tj0w SARS-CoV-2 and Influenza co-infection turned out to be a huge threat in recent times. The clinical presentation and disease severity is common in both the infection condition. The present paper deals with studying co-infection model system through systems biology approaches. Understanding signaling regulation in COVID-19 and co-infection model systems aid in the development of network-based models thereby suggesting intervention points for therapeutics. This paper highlights the aim of revealing such perturbations to decipher opportune mediating cross talks characterizing the deadly viral disease. The comparative analysis of both the models reveals major signaling protein NFκB and STAT1 playing a crucial role in establishing co-infection. By targeting these proteins at cellular level, it might help modulating the release of potent pro-inflammatory cytokines thereby taming the severity of the disease symptoms. Mathematical models developed here are precisely tailored and serves as a first step towards co-infection model offering flexibility and pitching towards therapeutic investigation. In 2019, coronavirus disease (COVID-19) originated in Wuhan, China and has been declared as global public health emergency [1, 2] . It is caused by Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) and the clinical symptoms includes dry cough, fever, shortness of breath, pharyngitis and tiredness [3] . When virus enters the lower respiratory tract, it shoots up the cytokine storm including IL6, IL12, IL1b, TNFa leading to viral sepsis, Acute Respiratory Distress Syndrome (ARDS) and multiple organ failure that ultimately leads to the death of the patient. The clinical symptoms and disease transmission mechanism of SARS-CoV-2 infection are similar to influenza virus infection [4] . During the early stage of the pandemic, one of the studies on 5700 patients hospitalized with COVID-19 infection in the New York City revealed that nearly 2.1% patients were found to have co-infection with other respiratory viruses including influenza A virus [5] . Furthermore, the percentage was increased in Wuhan, China where 4.35% confirmed cases of SARS-CoV-2 and influenza virus co-infection were reported [6] . Many parts of Europe and USA, considered it as 'Double threat' as its been anticipated in October 2020 that SARS-CoV-2 and Influenza virus might surge simultaneously [7] . One of the characteristic aspects of COVID-19 infection is the generation of cytokine storm. The molecular structure of SARS-CoV-2 has four major proteins namely SARS-CoV-2 spike (S), envelope (E), membrane (M) and nucleocapsid (N) [8] . Through its spike (S) glycoprotein, it interacts with alveolar epithelial type II cells which forms a protective layer on the inner respiratory tract. The high affinity interaction of Spike glycoprotein occurs with angiotensin-converting enzyme 2 (ACE2) and additionally processed by transmembrane serine protease 2 (TMPRSS2) both being present on the host cell [9] . The virus gets entry inside the cell, leading to the excessive production of pro-inflammatory cytokines such as IL6, IL1b, TNFa, IL12 and IL18. Most of these cytokines functions in acute inflammatory responses and are liable for causing Cytokine release syndrome (CRS) or Cytokine Storm syndrome [10] [11] [12] [13] . CRS is defined as cytokine-mediated systemic inflammatory response caused by excessive production of pro-inflammatory cytokines which ultimately results in clinical symptoms such as unrelenting high fever, lymphadenopathy, cytopaenia, central nervous system (CNS) abnormalities and if not treated leads to multiple organ failure (MOF) and death [14] . The hallmark of CRS is the unchecked feed forward activation and amplification of host cellular immune signaling which in turn activates the nearby cells to produce similar pattern of cytokines leading to cytokine storm [15] . Similar to COVID-19 infection, the cytokine storm plays a crucial role in severe influenza infected patients where aggressive pro-inflammatory response and in sufficient anti-inflammatory response leads to MOF and death of the patient [16] [17] [18] . Literature survey suggested that risk of death with SARS-CoV-2 was double in people which were prior infected with flu [19] . Furthermore, the current data suggested that the coinfection brings out synergistic effect of these two viruses which drastically increase the mortality rate as well as the demands of the health services [20] .This necessitated the design of new therapeutics targeting co-infection dynamics at basic cellular level and apparently help reduce the mortality rate. Macrophage plays crucial role in disease manifestation during both the infection. The activated alveolar epithelial cells release excessive amount of MCP-1 and GM-CSF which in turn recruits macrophages at the site of infection. The macrophages engulf the infected epithelial cells for phagocytosis. The single-cell transcriptomic profiling study of broncho alveolar lavage fluid samples from 88 patients with SARS-CoV-2-induced respiratory failure suggested that SARS-CoV-2 infects alveolar macrophages, which in turn respond by producing T cell chemo attractants [21] . Moreover, SARS-CoV-2 transcriptome has been found in alveolar macrophages. This shows that the capacity of SARS-CoV-2 to infect cells, is not only limited to alveolar epithelial cells but it can also infect alveolar macrophage. Further, macrophages play a central role in antiviral responses, tissue repair and fibrosis. Macrophages can be reprogrammed by environmental cues and thus can change their phenotype during an antiviral immune response as the viral infection progresses [22] . Thus, macrophages are considered as prominent candidate for therapeutics development in COVID-19 infections. The nature of both the viruses is similar. SARS-CoV-2 is positive strand, non-segmented RNA virus having their genetic material protected in viral envelope. The outermost layer has a spike glycoprotein which is recognised by the host cell whereas influenza is negative strand, segmented RNA virus, protected by viral envelope. The surface protein includes Haemagglutinin and Neuraminidase [23, 24] . Literature survey suggested that the RNA genome of SARS-CoV-1 and H1N1/H1N5 is recognized by Toll-like receptor 3 which is one of the major pattern recognition receptor (PRR) functioning in innate immune system [25] [26] [27] . Like TLR7, TLR8, TLR9 and TLR3 is present intracellular in an endosomal compartment and is known to recognize viral genome [28] . Consequently, TLR3 is likely to be involved in pathogenesis of SARS-CoV-2 virus infection. Further in SARS-CoV-1 infection, the stimulation of TLR3 activates NFjB transcription factor through TRIF signaling resulting in activation of pro-inflammatory genes such as IL6, TNFa, IL12p40 and IL1b [27] . Influenza A utilizes PI3K/ Akt pathway to activate transcription factor NFjB [29] . These cytokines are likely to activate nearby macrophages and epithelial cells resulting in massive production of cytokines and activation of macrophages, the phenomena called as Macrophage activated syndrome (MAS). One of the major genes activated by proinflammatory cytokines is Intercellular Adhesion Molecule (ICAM) and Vascular Cell Adhesion Molecule (VCAM). ICAM and VCAM function in cell-cell interaction. It is expressed by endothelial cell, epithelial cell and immune cells including macrophages. ICAM-1 interacts with LFA-1 present on leukocytes and facilitates its entry in tissues. ICAM-1 is found to enhance influenza viruses infection and their survival during an early stage of infection [30] . Similarly, VCAM also functions in leukocyte-endothelial cell interaction and both the cell adhesion molecules are likely to play a role in SARS-CoV-2 and influenza virus infection [31] [32] [33] . Together with IL6, TNFa, IL12 and IFNc, Macrophage Colony Stimulating Factor (MCSF) levels were also found to be elevated in COVID-19 infected patients [34] . MCSF is responsible for proliferation and differentiation of monocytes. MCSF plays a major role in establishing MAS during SARS-CoV-2 infection. Literature survey suggested the presence of curated population model defining spread of COVID-19 infection (e.g. BIOMD0000000957 or BIOMD0000000958) but is unable to define mathematically the co-infection dynamics. In the present work, using numerical simulation [35, 36] , combination of bits and pieces of literature evidences has acquired us to build new mathematical model which define the co-infection dynamics at cellular level and reveal crucial signaling proteins that can be targeted for generating new therapeutic regime for COVID-19 infection (Fig. 1) . For identifying each component in the interlinked signaling cascade, we have employed various immune signaling based databases. For IL6 signaling pathways we have used String database. TLR signaling pathway (map04620), PI3K-Akt signaling pathway (map04151), NFjB signaling pathway (map04064), TNFa signaling pathways (map04668), influenza A pathway (map05164) and SARS-CoV-2 infection pathway (map05171) has been obtained from KEGG pathway database [37] . For identifying targets of transcription factor NFjB, STAT1, STAT4 and STAT3 we have used TRRUST (version 2) database [38] . Here, we have constructed two mathematical models namely COVID-19 model and co-infection model (SARS-CoV-2 and influenza virus) using SimBiology toolbox from Matlab v7.11.1.866. SimBiology is a programmatic tool to model, simulate and analyse dynamic systems focusing on pharmacokinetic/pharmacodynam ics (PK/PD) and systems biology applications. It provides a block diagram editor for building models or one can create models pro-grammatically using the MATLAB Ò language. It supports systems biology mark-up language (SBML) format. Following are the steps used to create and simulate mathematical model of COVID-19 and co-infection model: The models were constructed in diagram editor using compartment, species, reactions and plot building blocks. In ''Model Building: Reaction Properties" we defined reactions and assigned specific rate laws (e.g. for association/ dissociation, we used Law of Mass action, for enzyme kinetics, we used Michaelis-Menten equation and for gene expression, we used Hill equation). Following this, concentration and parameter values are provided with respective units. As we are dealing with immune signaling, the concentration needs to be set between 10 3 -10 6 signaling molecules [39] . Further, the concentration of IL6 signaling pathway has been obtained from FACS experimental data [40] and has been incorporated to mimic the mathematical model close to physiological condition. The reactions were defined along with concentration and parameter followed by model simulation using ode15s (stiff/NDF) solver. The solver computes the model's state at the next time step using variable-order numerical differentiation formulas (NDFs). The values of absolute and relative tolerance are kept as default. The model is then simulated for 100 min close to 60 min of infection time taken by SARS-CoV-1 to infect the macrophages [41] . The behaviour of each component in the system have been obtained in the form of state v/s time graph. Parameter estimation is required to improve the accuracy of the mathematical model. The parameter of each reaction in reconstructed signaling network has been assigned within a range of parameter space. The simulation is performed to fine tune these parameters till a reproducible graph depicting the biochemical behaviour has been obtained. The process is iterative in nature as fine tuning of parameter space is required for proper estimation together with algorithmic data fitted to reach physiological accuracy. The mathematical model is analysed using various systems driven techniques such as Sensitivity analysis, Flux analysis, determining network properties and principal component analysis. The model is analysed through various aspects in order to churn out the crucial reactions that are governing the dynamics of the system (Fig. 2 ). Mathematically each biological network can be defined in terms of an interacting network of nodes and edges where biological components are the nodes and their interactions such as association, dissociation, phosphorylation, de-phosphorylation, ubiquitination, translocation etc. are the edges. Since the present work is based on immune signaling, the network is called as Protein-Protein Interaction (PPI) network. The network is analysed using Cytoscape (v 3.4.0) in terms of degree, clustering coefficient and centrality. The mathematical models were subjected to sensitivity analysis in order to check the robustness. Sensitivity analysis is a technique, wherein, the robustness of the mathematical model is achieved by perturbing the parameters of each reaction in the system. It further aids in parameter estimation close to physiological condition. Sensitivity analysis is obtained in terms of sensitivity coefficient and for the present work, we have performed local sensitivity analysis using SimBiology tool box from Matlab v7.11.1.866. We have calculated the time dependent sensitivities of the entire component in the system using SimBiology tool box. The SimBiology sensitivity analysis uses the ''complex-step approximation" to calculate derivatives of reaction rates. It calculates the time-dependent sensitivities of all the species states with respect to species initial conditions and parameter values in the model. Thus, if a model has species x and two parameters y, and z, then the time-dependent sensitivities of x with respect to each parameter value can be represented at the time-dependent derivatives as Eq. dx/dy, dx/dz where the numerators are the sensitivity output and the denominators are the sensitivity inputs to sensitivity analysis. Principal component analysis was performed to reduce the dimensionality of the large amount of sensitivity coefficients obtained from sensitivity analysis. The principal component (PC) score depicts the variation in sensitivities of each component in the system. The high PCA score depicts the high sensitivity of the component in the system or in other words the system might collapse if that component would be targeted. In network biology, highly connected nodes are those which has tendency to pass maximum biological information for one end to the other. The term ''collapse" here means, the transfer of biological information which is interrupted and which in turn greatly impacts the output of the system. PCA is meant for reducing dimensionality and removing the background noise of the biological data. For calculating PCA, we have used Matlab's function as score_cofficient = princomp A, where A = m*n matrix, sensitivity scores for each component with the other component in the system. The flux determines the productivity of each reaction, thereby determining the contribution of each reaction in the system. Comparative flux analysis is one of the optimum strategies to determine higher productive reactions in the system that may contribute to the disease pathogenesis. In the present work, we have obtained flux of each reaction by using COPASI pathway simulator (v4.11). The process is employed to reduce the complexity of the system by eliminating reaction with no or low impact on overall dynamics of the system. It simplifies the mathematical model through the elimination of redundant parameters which are not contributing to the system, thus easing in accurate prediction of the biological system. Here, we are employing sensitivity, flux and principal component analysis to churn out key components governing the dynamics of the system. The mathematical model has four compartments namely membrane, cytosol, nucleus and endosome. There are total 61 signaling components, 106 parameters, 58 reactions and kinetic laws (Fig. 3a) . The simulation was performed for 100 min, using Stiff Deterministic ODE15s solver (SimBiology toolbox) which generates the first order non-linear ODEs for each reaction in the system. The network shows production of four components at the end of simulation i.e. ICAM, VCAM, iNOS and MCSF (Fig. 3b) . The ODEs of the two major reactions are given below and for the entire model ODEs have been listed in supplementary file (S5). The ODEs represent the kinetic equations involving two major components of the system i.e. NFjB and STAT1, obtained after simulating the entire system for 100 min. The mathematical model has four compartments namely membrane, cytosol, nucleus and endosome. There are total 67 signaling components, 118 parameters, 64 reactions and kinetic laws (Fig. 4a) . The simulation was performed for 100 min, using Stiff Deterministic ODE15s solver (SimBiology toolbox) which generates the first order non-linear ODEs for each reaction in the system. The network shows production of four components at the end of simulation i.e. ICAM, VCAM, iNOS and MCSF (Fig. 4b) . The ODEs of two major reactions are given below (S6). The ODEs represents the kinetic equations involving two major components of the system i.e. NFjB and STAT1, obtained after simulating the co-infection By analysing both the models, we have deciphered the concentration of ICAM, VCAM, iNOS and MCSF, which is increased in coinfection model system as compared to COVID-19 model indicating the severity of the disease during co-infection (Fig. 5 ). The COVID-19 model consists of 66 nodes and 103 edges where as co-infection model consist of 72 nodes and 112 edges (Fig. 6) . The clustering coefficient of the COVID-19 network is 0.302 whereas co-infection network is 0.290, depicting that both the network are highly connected, therefore the flow of biological information is faster. Digging into the biology of the network, we have done network based ranking and identified five components which are found to have high degree as well as centrality value in both the network (Fig. 7) . Degree, Closeness, Betweenness centrality is the classic centrality measures. The betweenness centrality symbolizes the communication or flow of biological information in the system. In both the model systems, cytoplasmic and nuclear NFjB, phosphorylated STAT1, JAK, JAK1-JAK2 complex are the components through which maximum amount of biological information is transferred (Fig. 7a) . Closeness centrality determines the higher connectivity of the respective components in the system (Fig. 7b) and degree centrality determines the local connectivity of the components (Fig. 7c) . Cytoplasmic NFjB is found to be the cross talk identified between TNFa, IL1b, and TLR3 signaling pathway whereas phosphorylated STAT1 is the cross talk between IL6 and IFNc signaling pathway (Fig. 7d) . The flux of each reaction in both the models have been calculated in terms of molecules/ second and are compared to get broader picture of functioning of the two different models (Fig. 8a) . The comparative flux analysis shows that the reaction associated with major pro-inflammatory cytokines such as IL12, IL1b, IL6 and TNFa has high amount of flux (Fig. 8b & c) . This signifies that these reactions are responsible for major output of the system i.e., in the gene expression of ICAM, VCAM, MCP-1, iNOS etc. Targeting any one of the reactions may alter the expression levels of the above-mentioned gene resulting in decrease severity of the disease. Moreover, the flux of IL1b and IL6 reactions are increased which indicates their role in aggravated cytokine storm causing MAS (Fig. 8b & c) . The comparative analysis of PC score shows that major components common in both the systems are IL6, gp130, IL1b, VCAM, ICAM, NFjB and phosphorylated STAT1. The high score of IL6 and IL1b indicate their role in establishing pro-inflammatory immune response during infection stage. The data is supported by experimental data, where IL6 and IL1b levels has been found elevated in COVID-19 patients and IL6 is also considered as a parameter for measuring disease severity [42, 43] . Apart from this, PCA shows novel components such as cytoplasmic NFjB and phosphorylated STAT1. These components can be considered as targets as they have high centrality values. Targeting these components may affect the immune signaling established during SARS-CoV-2 and Co-infection. VCAM and ICAM is the end product of the system responsible for Macrophage Activated Syndrome (MAS). Principal component analysis shows the most sensitive or the most crucial components in both the model systems (Fig. 9) . The sensitivity coefficient of each component with respect to other components in the system have been documented in supplementary data (S7). The process of model reduction involves overlapping components and reactions identified in centrality, flux, sensitivity and principal component analysis. The complexities of both the models were reduced that resulted in 89.65% reduction in COVID-19 model whereas 90.62% reduction in co-infection model system. The reactions that are playing major role in regulating the dynamics of the system are as follows: The formation of cytokine storm is the hallmark for the severity of COVID-19 infection as well as influenza infection. The in-depth analysis of cellular signaling reveals several components that plays crucial role in disease progression. The high flux of IL6, IL1b, TNFa, IL12, IFNc signaling indicate the establishment of cytokine storm in COVID-19 and Co-infection model systems. After 100 min simulation, the end products are inflammatory proteins which include ICAM, VCAM, MCP-1, MCSF and iNOS. The co-infection model depicting the simultaneous infection of COVID-19 as well as influenza virus shows higher production of these inflammatory proteins indicating that the co-infection state is even worse than the COVID-19 infection as alone. The condition can be controlled once the productions of these inflammatory proteins are regulated at cellular level. Further, the model reduction reveals cytoplasmic NFjB and phosphorylated STAT1 as crucial therapeutic targets to regulate the intensity of cytokine storm established during infection. At present, there are various drugs available to inhibit the activity of NFjB such as Non-Steroidal Anti-inflammatory Drugs (NSAIDs). Aspirin and sodium salicylate are examples of antiinflammatory agents that target NFjB activity during chronic inflammation [44] . Cyclosporin A (CsA) and tacrolimus (FK-506) are immunosuppressive agents used during organ transplant and are found to inhibit calcineurin dependent NFjB activation [45] . The present line of therapeutics has their own potential side effects whereas use of peptide in therapeutical application provides more immunogenicity and specificity as well as has less off target effects. Its envisaged that the therapeutical design of the peptides is ideal to target the activity of NFjB for shorter duration of time. With the beginning of 2021, new strain of COVID-19 has emerged in United Kingdom [46] . The information on hostparasite interaction, proteins involved, mechanism of action is still unknown. It is foreseen that since we are targeting common cellu- lar proteins responsible for aggravated immune response the adopted strategy might be helpful or equally efficient against rising mutant strains by providing long-term immunity. The host immune signaling plays a crucial role in establishing or combating any viral infection. The present piece of work provides an in-depth system analysis of COVID-19 and Co-infection mathematical model which reveals major cellular protein (including transcription factors, cytokines, chemokines, receptors etc.) that can be targeted for future therapeutic intervention. The major components produced after 100 min simulations are IL6, IL12, IL1b, TNFa, IFNc, iNOS, VCAM, ICAM and MCSF. The production of iNOS, VCAM, ICAM and MCSF are high during Co-infection model indicating the disease severity during co-infection of influenza and SARS-CoV-2 virus. Further the network analysis highlighted the cellular protein JAK2, phosphorylated STAT1, cytoplasmic and nuclear NFjB showing high centrality value and indicating their role in establishing pro-inflammatory immune response during infection. These high centrality nodes are perfect targets to rewire the impaired immune response during SARS-CoV-2 infection. The models were further reduced to 90% through flux, sensitivity and principal component analysis and we have identified cytoplasmic proteins NFjB and phosphorylated STAT1as targets for future therapeutic intervention in case of both SARS-CoV-2 virus or Coinfection with the influenza virus. The underlying purpose here is to provide list of common potential targets out of numerous cellular proteins activated during SARS-CoV-2 infection and coinfection stage, which further might pave a way for generation of single therapeutics that is efficient in combating both the infection stages. The data generated here has opened new avenues for experimental investigations with accessible realms. Bhavnita Soni performed the mathematical modeling and was involved in acquisition, analysis and interpretation of the data. Shailza Singh was involved in the conception and design of study. Bhavnita Soni and Shailza Singh edited the manuscript and all authors approved the final version of the manuscript. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Clinical course and risk factors for mortality of adult inpatients with COVID-19 in Wuhan, China: a retrospective cohort study The Novel Coronavirus Originating in Wuhan, China Overview of immune response during SARS-CoV-2 infection: lessons from the past Comparing SARS-CoV-2 with SARS-CoV and influenza pandemics Presenting characteristics, comorbidities, and outcomes among 5700 patients hospitalized with COVID-19 in the The clinical characteristics of pneumonia patients coinfected with 2019 novel coronavirus and influenza virus in Wuhan Double threat of COVID-19 and influenza Structure, function, and antigenicity of the SARS-CoV-2 spike glycoprotein Unpuzzling COVID-19: tissue-related signaling pathways associated with SARS-CoV-2 infection and transmission The COVID-19 cytokine storm; what we know so far The pathogenesis and treatment of the 'Cytokine Storm' in COVID-19 COVID-19 cytokine storm: the anger of inflammation The cytokine storm in COVID-19: an overview of the involvement of the chemokine/chemokinereceptor system Cytokine storm syndrome in coronavirus disease 2019: a narrative review Cytokine storm and sepsis disease pathogenesis The cytokine storm of severe influenza and development of immunomodulatory therapy Outcomes of patients with severe influenza infection admitted to intensive care units: a retrospective study in a medical centre From influenza-induced acute lung injury to multiorgan failure Covid-19: Risk of death more than doubled in people who also had flu, English data show Interactions between SARS-CoV-2 and Influenza and the impact of coinfection on disease severity: a test negative design. medRxiv 2020 Circuits between infected macrophages and T cells in SARS-CoV-2 pneumonia Targeting macrophages as a therapeutic option in Coronavirus disease 2019 COVID-19 and flu, a perfect storm American society of microbiology. COVID-19 and the Flu Detrimental contribution of the toll-like receptor (TLR)3 to Influenza A virusinduced acute pneumonia Viral replication and innate host responses in primary human alveolar epithelial cells and alveolar macrophages infected with influenza H5N1 and H1N1 viruses Toll-like receptor 3 signaling via TRIF contributes to a protective innate immune response to severe acute respiratory syndrome coronavirus infection Effect of the phosphatidylinositol 3-kinase/Akt pathway on influenza A virus propagation ICAM-1 regulates the survival of influenza virus in lung epithelial cells during the early stages of infection Alveolar epithelial cells direct monocyte transepithelial migration upon influenza virus infection: impact of chemokines and adhesion molecules COVID-19, renin-angiotensin system and endothelial dysfunction Elevated expression of serum endothelial cell adhesion molecules in COVID-19 patients Clinical features of patients infected with 2019 novel coronavirus in Wuhan Stability analysis and numerical solutions of fractional order HIV/AIDS model Stability and numerical simulation of a fractional order plant-nectar-pollinator model KEGG: kyoto encyclopedia of genes and genomes TRRUST v2: an expanded reference database of human and mouse transcriptional regulatory interactions SnapShot: key numbers in biology Modeldriven experimental analysis of the function of SHP-2 in IL-6-induced Jak/STAT signaling Antibody-dependent infection of human macrophages by severe acute respiratory syndrome coronavirus Blocking IL-1 to prevent respiratory failure in COVID-19 The anti-inflammatory agents aspirin and salicylate inhibit the activity of IjB kinase-b A proteasome inhibitor prevents activation of NF-jB and stabilizes a newly phosphorylated form of IjB-a that is still bound to NF-jB About Variants of the Virus that Causes COVID-19 The authors would like to thank the Director, NCCS for supporting the Bioinformatics and High performance Computing Facility at National Centre for Cell Science (NCCS), PUNE, INDIA. Supplementary data to this article can be found online at https://doi.org/10.1016/j.csbj.2021.03.028.