Leishmaniasis is one of the major tropical and subtropical diseases caused by Leishmania spp which is a digenetic and dimorphic parasite residing inside the mammalian host. This disease is caused by female sand flies of Phlebotomus spp. and Lutzomyia spp. by bite and blood meal. Leishmaniasis is a disease that has been classified as an endemic throughout the world, affecting Asia, Africa, the Americas, and the Mediterranean region (Torres-Guerrero et al., 2017). In India, states such as Kerala, Madhya Pradesh, Haryana, Uttarakhand, Himachal Pradesh, Jammu, and Kashmir, Punjab, Assam, Delhi, Bihar, Rajasthan, Jharkhand, Uttar Pradesh and West Bengal has been majorly affected (Thakur et al., 2018) (Kumar et al., 2020). Over 90% of potentially fatal infections occur in these six countries viz., Brazil, Ethiopia, Sudan, South Sudan, India, and Bangladesh (Torres-Guerrero et al., 2017).

Broadly Leishmaniasis is classified into three forms viz., Cutaneous Leishmaniasis (CL), caused by Old World Leishmania spp such as Leishmania major/L. tropica and is characterized by the appearance of a variety of skin lesions, many of which are innocuous and self-healing. Mucocutaneous Leishmaniasis (MCL) is a long-term disease caused by New World Leishmania spp., such as L. braziliensis, that can cause extensive facial disfigurement and tissue destruction in the mouth and nose. The vast majority of MCL cases are not fatal. Complications associated with secondary infections, on the other hand, may result in death. Visceral Leishmaniasis (VL) infection caused by L. donovani/L. infantum is far more severe. It is frequently fatal if left untreated. It is caused by a systemic and progressive infection of macrophages in the reticulo-endothelial systems or lymphoid organs, primarily in the spleen, liver and bone marrow (Karunaweera and Ferreira, 2018).

Leishmania parasite has a digenetic life cycle completing in two stages. First, in the Phlebotomus sandfly vector and second in Mammalian host. When a sand fly takes a blood meal, it injects metacyclic promastigotes at the bite site. This stage of parasite is the infective stage with morphology characterized by slender shape and is elongated, motile and present extracellularly. Neutrophils are the first cells to be recruited at the bite site. Neutrophils engulf these promastigotes, but due to their short lifespan, they undergo apoptosis and are thus proposed to act as “Trojan horses” used by parasites to hitchhike ride in order to silently enter macrophages, avoiding cell activation (Regli et al., 2017). Inside the macrophage they divide and transform into amastigotes which infect the neighboring tissues. These amastigotes are small, circular and non-motile. When the sand fly takes blood meal from the infected individual, the infected macrophages are also taken up which burst and release amastigotes in the midgut of the sand fly. Later they mature and migrate to proboscis where they can infect another individual (Leta et al., 2014).

One of the prime inflammatory responses after infection is by our innate immunity. The parasite interacts with TLR2 which dimerizes with TLR6 and activates the MyD88 pathway. Through this pathway NFκB activates transcription of pro-inflammatory cytokines such as IL12, IFNγ, TNFα (Soni et al., 2016). These cytokines promotes the synthesis of Nitric oxide (NO) by iNOS which is directly involved in eliminating parasite by killing intracellular parasites within the NO-producing cell as well as those in bystander cells via diffusion across cell membranes (Scott and Novais, 2016). The response of pro-inflammatory cytokines is majorly controlled by IL10 as it is anti-inflammatory cytokine and has been reported to suppress the parasite killing and lead to its survival and persistence (Maspi et al., 2016). Hence, it becomes important to study these cytokines in a holistic way to understand the fundamental pathophysiology process of the disease.

Systems biology is a branch of biology that studies the complex mechanisms which underpin biological systems by treating gene, protein, biochemical network, and physiological responses as integrated parts of a larger system (Mc Auley et al., 2015). Being a holistic approach, Systems biology has been used to study various diseases for instance, using publicly available datasets as well as independent patient samples (Sambarey et al., 2017), identified a signature of ten genes that can distinguish tuberculosis from latent tuberculosis and HIV. In order to study COVID-19 and Influenza co-infection (Soni and Singh, 2021) used Systems biology approach and suggested therapeutic intervention points (Du and Elemento, 2015). stated Systems biology approaches are assisting in the understanding of tumour progression mechanisms and the development of more effective cancer therapies. Systems Biology approach versatility has thus created an urge to use it in studying cytokine reciprocity in Leishmaniasis and to identify players controlling this reciprocity to get a point of intervention; whose modulation may help in restoring the balance and shifting the paradigm to the desired phenotype.

Interaction of LPG (surface antigen of parasite) with TLR2 causes dimerization of TLR2 with TLR6 and activates MyD88 pathway leading to production of IL12 mediated by NFκB (Mukherjee et al., 2015). IL12 is known to induce IFNγ signaling through the JAK-STAT4 mediated pathway (Hamza et al., 2010). IFNγ further initiates a signaling cascade mediated by JAK-STAT1 and upregulates IL12 production (de Groen et al., 2015). Many transcription factors come into play for IL12 production such as IRF1, ICSBP, C/EBP and NFAT5 (Ma et al., 2015). IL12p40 promoter has Nucleosome1 whose histone modification is important in order to produce IL12. Rel family proteins can bind at Nucleosome 1 and lead to production of IL12, among which, NFAT5 is also a transcription factor from the same family (Albrecht et al., 2004) (López-Rodríguez et al., 2001). It can not only upregulate IL12 but also induce iNOS and NLRP3 inflammasome production (Buxadé et al., 2012) (Cornut et al., 2020) NFAT5 is activated by Brx mediated p38 pathway that often functions when the osmolarity of the cell increases, a phenomena which happens during inflammation and also in case of Leishmaniasis (Kino et al., 2009) (Bogdan, 2020). NFAT5 can not only promote parasite killing but it can bind at Sp1 binding site at IL10 promoter and inhibit its transcription (Choi et al., 2016).

Among all the transcription factors, NFAT5 drew our attention as it’s a transcription factor from the Rel family (Albrecht et al., 2004) (López-Rodríguez et al., 2001). TLR sensitive genes are activated by NFAT5. It does so by increasing the production of pro-inflammatory cytokines including IL6 and TNFα by recruiting NFκB and c-fos to their promoter regions. NFAT5 promotes demethylation of H3K27me3 marks from promoter region of pro-inflammatory genes and facilitates recruitment of other transcription factors such as NFκB and enhancers to accentuate their transcription. Sp1 is known to play a central role in IL10 expression (Brightbill et al., 2000; Lunazzi et al., 2021). Its binding at IL10 promoter occurs when H3 phosphorylation on Ser10 at Nucleosome2 takes place (Zhang et al., 2006). It lead us to the assumption and understanding that NFAT5 may induce chromatin remodeling at Nucleosome2 at IL10 promoter that consist of Sp1 binding site so that accessibility of Sp1 to IL10 Nucleosome2 region is hindered and IL10 transcription is hampered (Choi et al., 2016).

Parasite is notorious in nature. It secretes unknown effector molecules which induces chromatin remodeling in the host DNA (Afrin et al., 2019). The cell starts producing IL10 through ERK pathway and upon IL10 signaling the production of HDAC3 and NFIL3 (known IL12 inhibitors) increases leading to suppression of pro-inflammatory cytokine secretion and parasite survival (Yang et al., 2007) (Ma et al., 2015). Parasite also induces SHP-1 which inhibits NFAT5 by dephosphorylating it and thereby hampering parasite clearance and upregulating IL10 production (Zhou et al., 2014).

Using System’s Biology approach, we showed that NFAT5 can be an important transcription factor which can control IL10 and IL12 reciprocity at spatio-temporal level. It may do so by inducing chromatin remodeling at IL12p40 promoter and binding at Nucleosome1 region to upregulate IL12 and inhibiting IL10 synthesis by binding at Sp1 binding site on IL10 promoter and therefore inhibiting anti-inflammatory response. Thus, its presence/absence can play a vital role in deciding the fate of the disease. We also present NFAT5 as an essential component that can lead to parasite clearance by controlling the other determining parasite clearance machinery.

Common components in both the pathways after PCA were LPG, TLR2, TLR6, LPG-TLR2 and TLR6 complex, MyD88, IRAK4, IRAK1, TRAF6, IKK activation complex, NFκB, IL12, JAK2, 2STAT-1, NFAT5, iNOS, and IL10 (Figure 5). Other than the common components there were many other molecules which were model specific but also crucial as PCA allows identification and increase in the interpretability of critical components by minimizing information loss. Unique components which stood out from comparative analysis of both the model were NFIL3 and HDAC3 as they were specific to disease state model. IL10 can target NFIL3 to bind 10 kb upstream of the IL12p40 transcriptional start site and inhibit its synthesis. HDAC3 is involved in IL-10-induced histone deacetylation of the IL12p40 promoter (Ma X et al., 2015). The list of all the Principal components with their score is mentioned in Supplementary Material S2, S6.

Flux Analysis: As flux majorly contributes to the model phenotype; we observed following reactions had high flux:

According to the comparative flux analysis, the reaction associated with Healthy state model which have high flux are cytokines producing reactions such as IL12, IFNγ, osmolarity inducing reactions and iNOS inducing reaction. The reactions associated with diseased state which have high flux are inhibitory products of NFAT5 and IL12 suggestive of suppressing their synthesis. Hence, these high flux reactions can be associated with polarizing the model to healthy and diseased phenotype.

Model Reduction: For both the models we can observe the Quasi-potential landscapes (Figure 6). Both the graphs show us a dome shaped pattern or a peak like plot which shows the distribution of the model in the form of high PC, flux and components present at a higher concentration to fall at the top of the dome. As we move down the dome and peak; low PC, flux and concentration components will be present. As per PCA and flux; the reactions which contribute towards network output are the phenotype determinants.

Changes in the cytokines level in macrophages decide the fate of the disease. Reciprocity of cytokines have been previously investigated for IFN-γ and IL4 where, IFN-γ regulated Th2 response and IL4 regulated Th1 response during Leishmaniasis in murine models at temporal levels (KEMP et al., 1993). In murine leishmaniasis, IFN- γ and IL-4 are produced reciprocally in resolving or progressive infection, supporting the hypothesis that distinct T-helper subsets moderate the spectrum of this infectious disease (Heinzel et al., 1989). Moreover, L. donovani uses specific histone lysine methyltransferases/demethylases to redirect epigenetic programming of M(LPS + IFN-γ)/M(IL-10) genes for successful host establishment (Parmar et al., 2020). Hence, epigenetics have created an arena for regulating reciprocity between different cytokines, which makes cytokine reciprocity an enticing domain to study the remodeling of the immune response.

The results suggest that as the levels of IL12 decreases and IL10 increases the phenotype of the macrophage changes and that hampers the NO secretion which correspondingly inhibits parasite clearance and infectivity. Through our mathematical models we observed that when IL12 levels are high during the initial response, IL10 levels were low which promoted iNOS secretion and parasite clearance but as the time progressed the parasite promotes changes in the macrophages which hampers IL12 secretion and promotes IL10 secretion. As the secretion of IL10 increases, the levels of IL12, NFAT5 and iNOS decreases suggesting that parasite survival machinery has been activated (Figures 1-3). As, epigenetic regulation of host chromatin structure is influenced by the parasite (Afrin et al., 2019), our data shows that NFAT5 can be a key regulator of epigenetic modification of cytokines and may play a role in parasite elimination. NFAT5 presumably promotes inhibition of IL10 secretion by inhibiting H3Ser10 phosphorylation that recruits Sp1 on Nucelosome2 at the promoter region and at the same time it promotes demethylation of H3K27me3 marks from promoter region of IL12 at Nucleosome1 and recruits other transcription factors on IL12 promoter axis. These epigenetic modifications can be one of the most important changes responsible for changes in the immune response during Leishmaniasis (Krayem and Lipoldová, 2021). Over the course of time parasite that have skipped the clearance tends to hamper the ability of important transcription factors such as NFAT5, 2STAT-1 to participate in IL12 secretion by regulating them. Leishmania parasite is known to activate SHP-1 to dephosphorylate NFAT5 and prevent IL12 secretion. Moreover, parasite for its own survival modulates the host epigenetics by promoting ERK mediated H3 phosphorylation that makes Sp1 binding site again accessible to Sp1 for binding at IL10 promoter region (Yang et al., 2007) (Zheng et al., 2020). Since NFAT5 can now not bind at IL10 promoter by getting dephosphorylated by SHP-1, Sp1 can now access the modified H3 and carry out IL10 transcription (Figure 8). IL10 signaling promotes upregulation of IL12 inhibitors such as HDAC3 and NFIL3 which can now epigenetically and transcriptionally inactivate IL12 synthesis (Figures 9, 10).

Through PCA we tried to figure out key components which help in deciding the phenotype of the infected cell. The components of TLR2 and TLR6 mediated MyD88-NFκB pathway plays an important role in the phenotype decision making as they are known for promoting cytokine signaling (de Veer et al., 2003). Here, we also observed that IL12 and IL10 are principal components and they are of opposing nature hence, their synthesis and secretion levels will directly affect the parasite survival (Kane and Mosser, 2001) (Park et al., 2002). Jak2 is a tyrosine kinase which is recruited to IL12Rβ2, IFNγR2 to help in IL12 and IFNγ signaling (Bacon et al., 1995). iNOS is also a principal component as it is directly participating in parasite clearance by converting free N2 from metabolic pathways to NO (Olekhnovitch and Bousso, 2015). Another enticing aspect which came into the spotlight was transcription factors 2STAT-1 and NFAT5. Dimer form of phosphorylated STAT1 is well known in the literature for its ability to promote transcription of IFNγ (Bhardwaj et al., 2010). It is also known to activate Th1 response against parasites (Kima and Soong, 2013). NFAT5 is a member of the NFAT family of transcription factors which includes NFAT1, NFAT2, NFAT3 and NFAT4. These transcription factors are involved in osmoregulation of cells (Cen et al., 2020). It has been reported that NFAT5 binds to IL12p40 promoter region and accentuates IL12p40 synthesis in response to parasites (Buxadé et al., 2012; Tellechea et al., 2018). Since NFAT5 is also known to bind to the IL10 promoter region where Sp1 binding site is present and hampers IL10 synthesis (Choi et al., 2016), its dual functionality can not only help in increasing pro-inflammatory cytokine response but may also inhibit IL10 secretion so that parasite clearance can take place.

Flux Analysis decides the flow of molecules at a given time (Beguerisse-Díaz et al., 2018). In the healthy state, we observed that the reactions which were promoting IL12 secretion (Table 1) and the reactions that promoted IL10 production had low flux (S3 & S7). In the diseased state, the reactions which promoted IL10 synthesis and IL12 inhibition had high flux (Table 2). Thus the flux and PCA together lead us to figure out important components in the model which are (IL12, IL10, NFAT5, iNOS).

After performing model reduction the quasi potential landscape of both the models shows a dome and a peak like architecture. The reactions which were of high flux, high concentration and had high sensitivity coefficient lied at the top of the dome and peak whereas, the reactions with lower flux, concentration and sensitivity were at the base of the dome and peak (Figure 6). This materialized the crucial components and demonstrated that the model output or phenotype contributors could be IL12 and IL10 synthesis and secretion reactions, as well as reactions that inhibited their synthesis.

Crosstalk point helps in identification of the components which can connect two or more pathways. IL12 helps in IFNγ synthesis and its secretion is important to maintain healthy phenotype of the cells (Park et al., 2002), suggestive of the fact that for the parasite clearance IL12 is important. 2STAT-1 promotes synthesis of IFNγ as well as iNOS hence, STAT-1 is also promoting IL12 synthesis and promoting parasite clearance through iNOS upregulation (Bhardwaj et al., 2010) (Isnard et al., 2012) (Jayakumar et al., 2008). NFAT5 is not only promoting IL12 synthesis but also inhibiting IL10 (Choi et al., 2016). It also upregulates other processes which promote parasite clearance that is NLRP3 inflammasome production and iNOS secretion (Buxadé et al., 2012) (Cornut et al., 2020). Hence, NFAT5 has roles in promoting and inhibiting different processes which promote parasite survival. SHP-1 is a protein phosphatase which dephosphorylates NFAT5 (Yang et al., 2007) and therefore regulates NFAT5 activity in response to Leishmaniasis.

Through our computational analysis we could infer that NFAT5 may promote epigenetic remodeling in IL12 and IL10 promoter axis where it upregulates IL12 production and decreases IL10 secretion. Simultaneously, over the course of time as the parasite struggles for its survival, it upregulates SHP-1 in a manner that NFAT5 is dephosphorylated and inactivated for loss of its epigenetic modification regulating abilities, thereby changing the cytokines level altogether and apparently parasite survives to infect the other cells. Thus, NFAT5 can be taken as a point of intervention to altogether investigate the parasite clearance regime.

To conclude we posit that computational modeling of IL12 and IL10 reciprocity through NFAT5 in Leishmaniasis helps us to understand the signaling mechanism responsible for parasite survival. Thus, the system’s biological approach provides us with the necessary lead in a specific reconstructed biological network and by modulating NFAT5; we might be able to find a better therapeutic for treating Leishmaniasis.