Effects of Vitamin D Supplementation and Degradation on the Innate Immune System Response: Insights on SARS-CoV-2

Abstract

Vitamin D has been proven to be a strong stimulator of mechanisms associated with the elimination of pathogens. Because of its recognized effectiveness against viral infections, during SARS-CoV-2 infection, the effects of Vitamin D supplementation have been the object of debate. This study aims to contribute to this debate by the means of a qualitative phenomenological mathematical model in which the role of Vitamin D and its interactions with the innate immune system are explicitly considered. We show that Vitamin D influx and degradation can be considered as possible control parameters for the disease evaluation and recovery. By varying Vitamin D influx, three dynamical scenarios have been found with different modalities of recovery from the disease. Inside each scenario, Vitamin D degradation has been related to different degrees of severity in disease development. Interestingly, the emergence of hysteretic phenomenologies when Vitamin D influx is too low can be related to the onset of Long-COVID syndrome, confirming clinical evidence from recent studies on the topic.

1. Introduction

In 2019, mankind had to face the emergence of SARS-CoV-2, a new pathogen initially identified in Wuhan. This pathogen rapidly spread throughout the world, resulting in the COVID-19 pandemic and triggering an ongoing public health crisis. Confronting this pandemic presented a challenge that required massive material and intellectual resources. Today, after more than three years, we possess an enormous repository of scientific literature that attests to the efforts made to understand the virus and its mechanisms of action, with the hope of detecting effective ways to stop and control it.

COVID-19 displays a considerable diversity in its manifestation, ranging from predominantly asymptomatic and mild progressions to more severe and occasionally critical cases occurring in 10–20% of symptomatic individuals. This subset of patients faces a significant risk of mortality, exhibiting a range of affected organs and a variable array of symptoms [1]. The virus itself, the environment, the host, and their reciprocal interactions can be responsible for such heterogeneity [2].

To effectively react against disease, both weapons from the immune system are required: the innate immune system (including, among others, granulocytes, monocytes, neutrophils, macrophages) and the adaptive immune system (with T and B cells), see, i.e., [3]. The innate immune system, acting as the first line of defense, must be assigned a very special role in the protection against SARS-CoV-2 [4,5]. The innate immune system comprises physical defenses such as the epithelium, as well as chemical barriers including saliva and gastric acid. It also involves the complement cascade, antimicrobial peptides, and antigen-presenting cells such as dendritic cells and macrophages. In contrast to adaptive immunity, the responses of the innate immune system lack specificity and rely on a collection of phagocytic cells and proteins that identify shared characteristics of microorganisms. This rapid recognition leads to the swift elimination of infectious agents, encompassing viruses, bacteria, fungi, and parasites [6,7].

Recent findings have revealed that the innate immune system plays an important role in determining COVID-19 severity and outcome [8,9]. An insufficient innate immune response can trigger an elevated SARS-CoV-2 viral load, leading to an excessive activation of the adaptive immune system. This, in turn, prompts an intense production of pro-inflammatory cytokines, commonly known as the “cytokine storm,” which is an increase in C-reactive protein and D-dimer levels [10,11].

Although it is mainly recognized for its classical effects on musculoskeletal health, studies on nonclassical and extraskeletal actions of Vitamin D have increased exponentially [12,13]. In particular, during SARS-CoV-2 infection, many studies have shed new light on the role that Vitamin D may play in the fight against this virus. Several randomized clinical trials have demonstrated a significant association between its insufficiency and the severity as well as mortality of COVID-19. Moreover, these trials have suggested that supplementing it could potentially mitigate the disease’s seriousness [14]. Moreover, Vitamin D has been shown to regulate human immune functions by stimulating the innate immune response [15,16]. To delve into specifics, Vitamin D has been observed to bolster the protective physical barriers by reinforcing gap junctions, tight junctions, and adherens junctions that are mediated through E-cadherin. Moreover, it has been revealed to amplify the innate immune response by prompting the generation of numerous antimicrobial peptides such as defensins and cathelicidins within human monocytes and neutrophils. These peptides possess antibacterial, antiviral, and antifungal properties. Vitamin D also exerts an influence on cellular immunity. As a consequence of enhancing the innate immune system, Vitamin D can potentially lower the viral load and mitigate the excessive activation of the adaptive immune system. This, in turn, leads to the suppression of cytokine and C-reactive protein production, ultimately contributing to the reduction in COVID-19 severity and mortality; see, i.e., [17] and references therein.

This paper aims to delve into the intricate interplay between the innate immune system and Vitamin D in fighting pathogens. The immune system is one of the most complex, adaptive, and fascinating biological systems and many mathematical models have been created to explore its dynamics through many different approaches [18]. Among the most used are the theory of dynamical systems, where conventional integer-order differential equations are used to basically simulate the mechanistic functioning of the immune system [19,20,21]; the fractional-order modeling that turns out to be particularly appropriate to describe memory influences and hereditary properties of the immune systems [22,23,24]; and the game theory, where the pathogen and the immune system are considered as “players” that act systemically, by executing a series of strategies in a game situation that takes place on several rounds [25,26,27]. The present study moves exactly within the dynamical systems framework. On this line, drawing inspiration from a low-dimensional differential phenomenological model introduced by Young et al. [28], we explicitly consider the role of Vitamin D and its interactions with the innate immune system with the aim to enhance the understanding of how Vitamin D can influence the innate immune response against SARS-CoV-2.

In particular, the objective of this study is twofold: (i) to rationalize—by the means of a mathematical model—the role of Vitamin D in shaping the innate immune response against SARS-CoV-2, and (ii) to explore how this interaction can influence the severity and outcomes of COVID-19. Armed with this understanding, healthcare practitioners could better navigate the complex landscape of the disease and formulate more targeted strategies for intervention. This would aid not only in the direction of refining COVID-19 treatment therapeutic modalities but would also align the larger goal to prevent its severe manifestations, potentially alleviating the burden on healthcare systems and improving patient well-being.

The novelty of this research lies in its innovative approach to investigate the role of Vitamin D in the context of COVID-19. The integration of the innate immune system and Vitamin D into a comprehensive mathematical model is in fact a feature that enriches the existing literature on the subject since (to the best of our knowledge) no other mathematical model has been used to rationalize and investigate the mechanisms underlying their reciprocal interplay, also demonstrating the possibility of different outcomes as a result of this specific interaction.

This paper is structured as follows. Section 2 presents experimental findings derived from recent studies that emphasize the impacts of Vitamin D in relation to the treatment of COVID-19. These findings serve as inspiration for our current research and are utilized to validate the analytical outcomes we achieve. In Section 3, we introduce the model, detect its equilibria, their stability properties, and the bifurcations involved. The possible emergence of a hysteretic phenomenology because of the occurrence of a backward scenario has been investigated. In Section 4, the consequences of the backward and forward bifurcations, in terms of possible dynamical scenarios, are discussed. In Section 5, some conclusions end the paper.

2. Vitamin D and COVID-19: Experimental Evidences

The role of Vitamin D in the prevention and treatment of COVID-19 has recently been at the center of an intense scientific debate. This issue has become of even greater importance because of the potential negative impact of hypovitaminosis D on the incidence of SARS-CoV-2 infection and on the prognosis of COVID-19 [29]. However, how to treat and manage the Vitamin D and COVID-19 relationship is currently the object of investigations.

In this respect, the Glucocorticoid Induced Osteoporosis Skeletal Endocrinology Group commissioned a group of expert shareholders that, by a critical and systematic review of the recent literature on the subject, highlighted some fixed points in this vexata quaestio [30]. For example, deficient levels of Vitamin D associated with hypocalcemia were frequently described also in patients severely affected by other pandemics such as Ebola in 2016 and SARS in 2003 [31]; this finding was also present in about 80% of patients hospitalized in Italy for SARS-CoV-2 [32]. Indeed, the lower the serum calcium level, the worse the prognosis [33]. Moreover, in a medium-long period of observation, a low level of Vitamin D turned out to be decisive in influencing the patient’s comorbidities and might therefore, albeit in a variable way, be linked to the severity of the COVID-19 [34]. Furthermore, a noticeable correlation between disease outbreaks and geographic latitude was identified, indicating a connection between Vitamin D levels and reduced sun exposure [35,36]. This relationship becomes evident in the case of Italy, a European nation characterized by a considerable prevalence of hypovitaminosis D. Intriguingly, Italy also recorded the highest incidence of SARS-CoV-2 and COVID-19 infections, with the northern regions experiencing the most pronounced impact [37]. A correlation between Vitamin D deficiency and the severity of SARS-CoV-2 infections was observed above all in the elderly, who are often suffering from multiple chronic pathologies that can affect the severity of the infection. More precisely, 65% of the most severe cases also had a Vitamin D deficiency compared to the proportion of patients who contracted COVID-19 in a milder form. Vitamin D deficiency was also related to higher hospitalization and higher mortality from COVID-19 [38]. Furthermore, evidence also highlights the positive outcomes associated with supplementing Vitamin D during the progression of the infection. In a retrospective analysis, the use of cholecalciferol as a supplementary therapy was linked to a decreased risk of mortality from COVID-19 [39]. Moreover, the administration of a substantial dosage of calcifediol was observed to lower the requirement for care in the Intensive Care Unit (ICU) among patients admitted to the hospital. The multivariate assessment of risk for ICU patients who were or were not given calcifediol yielded a striking figure of $0.03$ (with a confidence interval of 95%, ranging from $0.003$ to $0.25$) [40]. Just recently, the favorable effects of Vitamin D supplementation in individuals hospitalized due to COVID-19 received further endorsement through a meticulously designed randomized, double-blind, placebo-controlled trial [41].

To further elucidate the role of Vitamin D in the pathogen control, mathematical modeling stands as a valuable tool due to its potential to offer potential scenarios associated with the implementation of particular clinical measures. In this respect, mathematical models can act at two different and somehow complementary levels. Those that take into account as many details as possible of the complex ’world’ related to the immune system functioning are usually described by sophisticated mathematical tools and can give a contribution in terms of quantitative predictions. In contrast, qualitative phenomenological models capitalize on a more straightforward mathematical approach and can provide a deep mechanistic comprehension of the underlying mechanisms. These models emphasize the roles and interconnections of these mechanisms in a profound manner.

Basing on this latter approach and with the aim to gain a deeper understanding of the empirical evidence discussed above, in the next section we introduce and investigate a low-dimensional phenomenological model in which the role of Vitamin D and its interactions with the innate immune system is explicitly included.

3. The Model

To rationalize the possible role of Vitamin D in the pathogen control, let us introduce the following low-dimensional model:

$$\begin{array}{ccc}{\displaystyle \frac{dx}{dt}}& =\hfill & rx-\beta xz,\hfill \\ \\ {\displaystyle \frac{dz}{dt}}& =\hfill & \delta x-\gamma xz+\eta -\mu z+kD,\hfill \\ \\ {\displaystyle \frac{dD}{dt}}& =\hfill & \mathsf{\Gamma }-hD-\omega xD,\hfill \end{array}$$

(1)

where $x\left(t\right)$ is the concentration of pathogen load; $z\left(t\right)$ is the concentration of active PolyMorphonuclear Neutrophils (PMNs), namely neutrophils and other rapidly responding elements of the innate immune system and D(t) represents the concentration of Vitamin D. In the equation for the pathogen dynamics, the term $rx$ describes the pathogen growth whereas $-\beta xz$ represents the phagocytic action of the PMNs. In the equation for the dynamics of the active PMNs, the term $\delta x$ represents the influx of PMNs into the alveolar sacs as dependent only on the viral load x whereas $-\gamma xz$ indicates that PMNs die after ingesting and destroying pathogens. The constant term $\eta$ is a small baseline influx of PMNs; $-\mu z$ is the natural death term of PMNs whereas $kD$ points out that Vitamin D spurs PMNs activation. In the equation for the Vitamin D dynamics, Vitamin D is assumed to be supplied at a constant rate $\mathsf{\Gamma }$, from dietary and external sources; $-hD$ represents the degradation of Vitamin D whereas the term $-\omega xD$ represents the degradation of Vitamin D due to the presence of pathogen. All the parameters in model (1) are assumed to be positive constants. For a detailed description of each parameter, please refer to

Table 1. Model parameters, description and values chosen for the discussions.

Parameter	Description	Value
$\gamma$	neutrophils apoptosis rate	$0.5$
$\eta$	neutrophils recruitment rate	1
$\mu$	neutrophils death rate	$1.2$
$\delta$	neutrophils elicitation rate	2
$\beta$	neutrophils phagocytosis rate	$0.9$
r	pathogen growth rate	$0.8$
$\omega$	Vitamin D degradation rate due to pathogen	$0.5$
k	neutrophils activation rate by Vitamin D	$0.3$
$\mathsf{\Gamma }$	Vitamin D supply rate	varying
h	Vitamin D degradation rate	varying

.

In the following, we carry on our investigations under the assumptions that: (i) $\beta >\frac{\gamma r}{\delta }$, meaning that a large pathogen load may not be able to overcome the innate system; (ii) $\frac{\eta }{\mu }\ll \frac{\delta }{\gamma }$, meaning that the background level of PMNs is much less than the effective carrying capacity.

In the next section, we will investigate on the positivity and boundness of solutions as well as on the existence of equilibria for model (1).

3.1. Positivity and Boundness of Solutions

Let us consider strictly positive initial conditions for model (1), i.e., $x\left(0\right)>0$, $z\left(0\right)>0$ and $D\left(0\right)>0$. By direct integration of the first equation of model (1), we observe that

$$x\left(t\right)=x\left(0\right)exp\left[{\int }_0^t\left(r-\beta z\left(s\right)\right)ds\right],$$

so that $x\left(t\right)$ solution has the same sign of the initial condition, i.e., it is strictly positive. Taking this into account, by the third equation in model (1), it follows $\frac{dD}{dt}\ge -D(h+\omega x)$, so that

$$D\left(t\right)\ge D\left(0\right)exp\left[-{\int }_0^t\left(h+\omega x\left(s\right)\right)ds\right],$$

and $D\left(t\right)$ is strictly positive. Finally, because of the positivity of $x\left(t\right)$ and $D\left(t\right)$, from the second equation in model (1), it follows $\frac{dz}{dt}\ge -z(\gamma x+\mu )$ so that

$$z\left(t\right)\ge z\left(0\right)exp\left[-{\int }_0^t\left(\mu +\gamma x\left(s\right)\right)ds\right],$$

and $z\left(t\right)$ is strictly positive. Therefore the non negative orthant is positively invariant for system (1).

As far as boundness is concerned, we recall that our assumptions on the model parameters implies the following inequalities to hold: $r/\beta <\delta /\gamma$ and $\eta /\mu \ll \delta /\gamma$. In addition, because of the positivity of the solutions, from the third equation in system (1) one has $\frac{dD}{dt}\le \mathsf{\Gamma }-hD,$ so that ${lim\; sup}_{t\to +\infty }D\left(t\right)\le \mathsf{\Gamma }/h$. Therefore the upper bound for Vitamin D is given by the ratio between its supply and its degradation.

We also observe that, in the absence of pathogen ($x=0$), it can be easily seen that the $z\left(t\right)$ and $D\left(t\right)$ variables in system (1) are bounded for all $t\ge 0$. When the pathogen is introduced ($x\ne 0$) and start increasing, then the neutrophil variable $z\left(t\right)$ increases too, hence limiting the exponential increasing of the pathogen $x\left(t\right)$.

Suppose that a $t_0>0$ exists such that $z\left(t\right)\ge r/\beta$ for $t\in [t_0,T]$, with $t_0<T\le \infty$. Therefore in $[t_0,T]$, from the first equation of model (1), we have that $\frac{dx}{dt}=x(r-\beta z)<0$ so that $x\left(t\right)$ is decreasing in $[t_0,T]$ and a $M>0$ exists such that $x\left(t\right)\le x\left(t_0\right)=M$, $\forall t\in [t_0,T]$. Now, let us consider the auxiliary variable $W\left(t\right)=x\left(t\right)+z\left(t\right)$, $\forall t\in [t_0,T]$. It follows that:

$${\displaystyle \frac{dW}{dt}}={\displaystyle \frac{dx}{dt}}+{\displaystyle \frac{dz}{dt}}=x(r-\beta z)+x(\delta -\gamma z)+\eta -\mu z+kD\le \eta +(\delta +\mu )x+kD-\mu (x+z),$$

so that ${\displaystyle \frac{dW}{dt}}+\mu W\le {\mathsf{\Sigma }}_{\ast },$ with ${\mathsf{\Sigma }}_{\ast }=\eta +(\mu +\delta )M+k\mathsf{\Gamma }/h$. This leads to:

$${\displaystyle W\left(t\right)\le W\left(t_0\right)e^{-\mu (t-t_0)}+\frac{{\mathsf{\Sigma }}_{\ast }}{\mu }\left(1-e^{-\mu (t-t_0)}\right).}$$

Therefore ${\displaystyle \underset{t\to +\infty }{lim\; sup}W\left(t\right)\le {\mathsf{\Sigma }}_{\ast }/\mu }$, so that $x\left(t\right)$ and $z\left(t\right)$ are bounded $\forall t\in [t_0,T]$. Such an upper bound is then given by the background level of $PMNs$ increased by a factor proportional to the pathogen level at the time $t_0$ and by a factor proportional to the ratio between the vitamin D supply and its degradation.

Remark 1. 

If we suppose $z\left(t\right)<r/\beta$ for $t\in [0,\tilde{t}]$, then by the first and second equation in model (1), it follows that both $\frac{dx}{dt}>0$ and $\frac{dz}{dt}>0$. Therefore a $t_0>\tilde{t}$ will exist such that $z\left(t\right)\ge r/\beta$ in $[t_0,T]$ so that we come back to the previous case.

3.2. Equilibria and Their Properties

System (1) admits the pathogen-free equilibrium $P_0=(x_e,z_e,D_e)$ where,

$${\displaystyle x_e=0,z_e=\frac{1}{\mu }\left({\displaystyle \eta +k\frac{\mathsf{\Gamma }}{h}}\right),D_e=\frac{\mathsf{\Gamma }}{h}.}$$

(2)

It may also admit pathogen-persistent equilibria $P^{\ast }=(x^{\ast },z^{\ast },D^{\ast })$ where

$${\displaystyle z^{\ast }=\frac{r}{\beta },D^{\ast }=\frac{\mathsf{\Gamma }}{h+\omega x^{\ast }},}$$

(3)

with $x^{\ast }$ positive root of the second order algebraic equation

$$A{x}^2+Bx+C=0.$$

(4)

Here,

$$\begin{array}{c}A=\omega \left(\beta \delta -\gamma r\right),\hfill \\ \\ B=h\left(\beta \delta -\gamma r\right)+\omega \left(\beta \eta -\mu r\right),\hfill \\ \\ C=\mathsf{\Gamma }\beta k+h\left(\beta \eta -\mu r\right).\hfill \end{array}$$

(5)

Because of the assumption ${\displaystyle \beta >\frac{\gamma r}{\delta }}$, it follows $A>0$. We also observe that if $\beta \eta >\mu r$, then $B>0$ and $C>0$. Therefore, no positive pathogen-persistent equilibrium can exist in this case.

Let us introduce the following quantities:

$${\displaystyle h_1=\frac{\omega (\mu r-\beta \eta )}{\beta \delta -\gamma r};h_2=\frac{\mathsf{\Gamma }\beta k}{\mu r-\beta \eta },}$$

(6)

and

$${\displaystyle {\mathsf{\Gamma }}_b=\frac{\omega {\left(\mu r-\beta \eta \right)}^2}{\beta k\left(\beta \delta -\gamma r\right)}.}$$

(7)

In the next, we also assume the following assumption:

$${\displaystyle \frac{\gamma r}{\delta }<\beta <\frac{\mu r}{\eta }.}$$

(8)

We are hence in the position to state the following results:

Proposition 1. 

Let inequalities (8) hold. Then:

(i) 

$B>0⟺h>h_1$.

(ii) 

$C>0⟺h<h_2$.

(iii) 

$h_1<h_2⟺\mathsf{\Gamma }>{\mathsf{\Gamma }}_b$.

Proof. 

Results (i)–(iii) easily follow by direct calculation. □

Proposition 2. 

Let inequalities (8) hold. Then, the second order algebraic Equation (4) admits real roots if and only if $h<h_2^{\ast }\vee h>h_1^{\ast }$, where

$${\displaystyle h_{1/2}^{\ast }=\frac{\omega \left(\mu r-\beta \eta \right)\mp 2\sqrt{{\mathsf{\Delta }}^{\ast }}}{\gamma r-\beta \delta },}$$

(9)

and

$${\mathsf{\Delta }}^{\ast }=\mathsf{\Gamma }\beta k\omega \left(\beta \delta -\gamma r\right).$$

(10)

Proof. 

The second-order algebraic Equation (4) for the variable x has real roots if and only if its discriminant $B^2-4AC>0$. Let’s define $a=\beta \delta -\gamma r$. Using the condition given by (8), we establish that $a>0$ and ${\mathsf{\Delta }}^{\ast }=\mathsf{\Gamma }\beta k\omega a>0$. Furthermore, if we set $b=\mu r-\beta \eta$, direct calculation yields:

$$\begin{array}{cc}\hfill B^2-4AC& =a^2{h}^2+2\omega abh+{\omega }^2{b}^2-4\mathsf{\Gamma }\beta k\omega a\hfill \\ \hfill & =\mathrm{a}\mathrm{second}-\mathrm{order}\mathrm{polynomial}\mathrm{in}h.\hfill \end{array}$$

Denoting the discriminant of this polynomial as $\mathsf{\Delta }$, we find that:

$$\begin{array}{cc}\hfill \sqrt{\frac{\mathsf{\Delta }}{4}}& =\sqrt{{\omega }^2{a}^2{b}^2-a^2({\omega }^2{b}^2-4\mathsf{\Gamma }\beta k\omega a)}=2a\sqrt{\mathsf{\Gamma }\beta k\omega a}=2a\sqrt{{\mathsf{\Delta }}^{\ast }},\hfill \end{array}$$

which implies that its roots are real and given by:

$$\begin{array}{cc}\hfill h_{1/2}^{\ast }& =\frac{-\omega ab\pm 2a\sqrt{{\mathsf{\Delta }}^{\ast }}}{a^2}=\frac{-\omega b\pm 2\sqrt{{\mathsf{\Delta }}^{\ast }}}{a}=\frac{\omega b\mp 2\sqrt{{\mathsf{\Delta }}^{\ast }}}{-a}.\hfill \end{array}$$

From (8), we deduce that $b>0$, ensuring that ${\displaystyle h_2^{\ast }=\frac{\omega b+2\sqrt{{\mathsf{\Delta }}^{\ast }}}{-a}<0}$, and $h_2^{\ast }<h_1^{\ast }$, where $h_1^{\ast }$ can be either positive or negative.

Finally, we establish:

$$B^2-4AC>0⟺h<h_2^{\ast }\vee h>h_1^{\ast }.$$

□

Moreover, the following proposition provides useful results on the above threshold quantities:

Proposition 3. 

Let inequalities (8) hold. Then:

(i) 

${\displaystyle h_1^{\ast }>0⟺\mathsf{\Gamma }>\frac{1}{4}{\mathsf{\Gamma }}_b}$.

(ii) 

$h_1>h_1^{\ast }⟺\mathsf{\Gamma }<{\mathsf{\Gamma }}_b$.

(iii) 

$h_2>h_1^{\ast }\forall \mathsf{\Gamma }\ne {\mathsf{\Gamma }}_b$.

Proof. 

Results (i)–(iii) follow by direct calculation. In fact:

(i) $h_1^{\ast }>0⟺2\sqrt{{\mathsf{\Delta }}^{\ast }}>\omega \left(\mu r-\beta \eta \right)⟺4\mathsf{\Gamma }\beta k\left(\beta \delta -\gamma r\right)>\omega {\left(\mu r-\beta \eta \right)}^2$.

Therefore $h_1^{\ast }>0⟺\mathsf{\Gamma }>\frac{1}{4}{\mathsf{\Gamma }}_b$.

(ii) $h_1>h_1^{\ast }⟺\mathsf{\Gamma }\beta k\left(\beta \delta -\gamma r\right)<\omega {\left(\mu r-\beta \eta \right)}^2⟺\mathsf{\Gamma }<{\mathsf{\Gamma }}_b$.

(iii) $h_2-h_1^{\ast }>0⟺\mathsf{\Gamma }\beta k(\beta \delta -\gamma r)+\omega {(\beta \eta -\mu r)}^2-2(\mu r-\beta \eta )\sqrt{\mathsf{\Gamma }\beta k\omega \left(\beta \delta -\gamma r\right)}$>0.

Therefore $h_2-h_1^{\ast }>0\forall \mathsf{\Gamma }\ne {\mathsf{\Gamma }}_b$. □

By combining the above results we can establish conditions, in terms of the system parameters, for the existence of the pathogen-free and pathogen persistent equilibria:

Theorem 1. 

Let inequalities (8) hold. Then:

(i) If ${\displaystyle 0<\mathsf{\Gamma }<\frac{1}{4}{\mathsf{\Gamma }}_b}$, model (1): admits a pathogen-free equilibrium $\forall h>0$; admits two pathogen-persistent equilibria $\forall h\in \left[0,h_2\right]$; admits one pathogen-persistent equilibrium for $h>h_2$.

(ii) If ${\displaystyle \frac{1}{4}{\mathsf{\Gamma }}_b<\mathsf{\Gamma }<{\mathsf{\Gamma }}_b}$, model (1): admits a pathogen-free equilibrium $\forall h>0$; does not admit pathogen-persistent equilibria $\forall h\in \left[0,h_1^{\ast }\right]$; admits two pathogen-persistent equilibria $\forall h\in \left[h_1^{\ast },h_2\right]$; admits one pathogen-persistent equilibrium for $h>h_2$.

(iii) If $\mathsf{\Gamma }>{\mathsf{\Gamma }}_b$, model (1): admits a pathogen-free equilibrium $\forall h>0$; does not admit pathogen-persistent equilibria $\forall h\in \left[0,h_2\right]$; admits one pathogen-persistent equilibrium for $h>h_2$.

Proof. 

It is easy to observe that model (1) admits the pathogen-free equilibrium $P_0$ for every value of the parameters. Results on the existence of the pathogen-persistent equilibria follow by the study of the roots of the second order algebraic Equation (4) by the Descartes’ rule of signs and by properly combining results given in Propositions 1–3. In fact:

(i) If ${\displaystyle 0<\mathsf{\Gamma }<\frac{1}{4}{\mathsf{\Gamma }}_b}$, then the following relation holds: $h_2^{\ast }<h_1^{\ast }<0<h_2<h_1$, so that two real positive roots exist for $0<h<h_2$ and one real positive root exists for $h>h_2$.

(ii) If ${\displaystyle \frac{1}{4}{\mathsf{\Gamma }}_b<\mathsf{\Gamma }<{\mathsf{\Gamma }}_b}$, then the following relation holds: $h_2^{\ast }<0<h_1^{\ast }<h_2<h_1$, so that no real root exists for $0<h<h_1^{\ast }$; two real positive roots exist for $h_1^{\ast }<h<h_2$ and one real positive root exists for $h>h_2$.

(iii) If $\mathsf{\Gamma }>{\mathsf{\Gamma }}_b$, then the following relation holds $h_2^{\ast }<0<h_1<h_1^{\ast }<h_2$, so that no real root exists for $0<h<h_1^{\ast }$; no real positive root exists for $h_1^{\ast }<h<h_2$ and one real positive root exists for $h>h_2$. □

These findings turn out to be compatible with the occurrence of a backward scenario. In the next subsection, we define the circumstances under which the model showcases either a forward or a backward phenomenology, and we assess the plausible influence of Vitamin D on both of these scenarios. Specifically, the backward phenomenology is tied to a state of bistability between the equilibrium where the pathogen is absent and the equilibrium where the pathogen persists. This condition could drive the system toward behaviors resembling hysteresis. Broadly speaking, the term “hysteresis” conveys the notion of “irreversibility,” indicating enduring effects that linger even after the triggering factors have been eliminated. Its link with hysteresis is the reason why many papers, even in different contexts, have been focused on backward bifurcation, i.e., [42,43,44,45].

3.3. Forward versus Backward Bifurcation Scenarios

In mathematical epidemiology, much attention has been specifically devoted to the backward bifurcation [46,47,48,49,50]. Within this context, the reference control quantity is the basic reproduction number $R_0$, a metric commonly characterized as the expected number of new infections produced by a single infective individual introduced into a disease-free population [51] Notably, when $R_0=1$, this signifies the critical threshold where the disease-free equilibrium transitions between states of stability and instability.

At $R_0=1$, two distinct bifurcation scenarios can arise: (i) a forward bifurcation emerges, signifying the possibility of disease elimination when $R_0$ drops below the threshold of $R_0=1$; (ii) a backward bifurcation scenario unfolds, characterized by a saddle-node (sn) bifurcation occurring at $R_0=R_0^{sn}<1$, accompanied by a subcritical transcritical bifurcation at $R_0=1$. This scenario introduces a multitude of endemic equilibria and a subcritical persistence of the disease. Consequently, eradicating the disease necessitates not only reducing $R_0$ below 1 but lowering it below the critical value $R_{sn}$. From the perspective of disease control, the identification and effective management of backward bifurcations hold significant importance [43,44,48,52,53,54].

To detect the occurrence of a backward phenomenology, we start by investigating the stability properties of the pathogen-free equilibrium $P_0$. In this respect, we recall that model (1) always admits the pathogen-free equilibrium $P_0$. As far as its stability properties are concerned, the following result holds:

Theorem 2. 

The pathogen-free equilibrium $P_0$ is locally asymptotically stable for $h<h_2$ and unstable for $h>h_2$.

Proof. 

The Jacobian matrix of the model (1), evaluated at point $P_0$, takes the form:

$$J\left(P_0\right)=\left[\begin{array}{ccc}{\displaystyle -\beta \frac{\mathsf{\Gamma }k+\eta h}{h\mu }+r}& 0& 0\\ {\displaystyle -\gamma \frac{\mathsf{\Gamma }k+\eta h}{h\mu }+\delta }& -\mu & k\\ {\displaystyle -\frac{\mathsf{\Gamma }\omega }{h}}& 0& -h\end{array}\right],$$

(11)

This matrix admits as eigenvalues ${\lambda }_1=-\mu <0$, ${\lambda }_1=-h<0$ and ${\displaystyle {\lambda }_3=-\beta \frac{\mathsf{\Gamma }k+\eta h}{h\mu }+r}$. It is straightforward to demonstrate that:

$${\lambda }_3<0⟺r<\beta \frac{\mathsf{\Gamma }k+\eta h}{h\mu }⟺h(\mu r-\beta \eta )<\mathsf{\Gamma }\beta k⟺h<\frac{\mathsf{\Gamma }\beta k}{\mu r-\beta \eta }.$$

Thus, point $P_0$ is locally asymptotically stable when $h<h_2$, where $h_2$ is defined in (6). Conversely, it becomes unstable when h exceeds this threshold. □

By Theorems 1 and 2, a backward-type scenario can be glimpsed when $0<\mathsf{\Gamma }<{\mathsf{\Gamma }}_b$, with $h=h_2$ as a transcritical threshold. As stated in the following Theorem, it is the level of Vitamin D supplementation, by the mean of the parameter $\mathsf{\Gamma }$, that can make the difference between the forward and the backward scenario:

Theorem 3. 

Let inequalities (8) hold. Then: (i) If $\mathsf{\Gamma }<{\mathsf{\Gamma }}_b$, then system (1) exhibits a backward bifurcation at $h=h_2$. (ii) If $\mathsf{\Gamma }>{\mathsf{\Gamma }}_b$ then system (1) exhibits a forward bifurcation at $h=h_2$.

Proof. 

The determination of whether the backward or forward scenario takes place is examined through the approach outlined in [42]. The coefficients within the equilibrium Equation (4) can be seen as being dependent on the parameter h. Consequently, at $h=h_2$, it holds that $C\left(h_2\right)=0$, and as a result, Equation (4) transforms into:

$$A\left(h_2\right)x^2+B\left(h_2\right)x=0,$$

with roots $x=0$ and $x=-\frac{B\left(h_2\right)}{A\left(h_2\right)}$. The root $x=0$ corresponds to the pathogen-free equilibrium, while $x=-\frac{B\left(h_2\right)}{A\left(h_2\right)}$ represents a positive pathogen-persistent equilibrium, but only under the condition that $B\left(h_2\right)$ and $A\left(h_2\right)$ exhibit opposite signs. Consequently, since it’s established that $A>0$, it follows that $B\left(h_2\right)<0$ must be satisfied.

Now, when Equation (4) is implicitly differentiated with respect to h, we obtain:

$${\displaystyle (2Ax+B)\frac{dx}{dh}+\frac{dA}{dh}{x}^2+\frac{dB}{dh}x+\frac{dC}{dh}=0,}$$

that evaluated at $x=0$ and $h=h_2$, returns:

$${\displaystyle B\left(h_2\right)\frac{dx}{dh}\left(h_2\right)=-\frac{dC}{dh}\left(h_2\right)>0,}$$

(12)

since, by (5), one has that ${\displaystyle \frac{dC}{dh}\left(h_2\right)=\beta \eta -\mu r<0}$. Therefore the slope of the bifurcation curve at $x=0$ must have the same sign of $B\left(h_2\right)$. As a consequence: (i) if $B\left(h_2\right)<0$ then a backward bifurcation occurs at $h=h_2$; (ii) if $B\left(h_2\right)>0$, then the system displays a forward bifurcation at $h=h_2$. By (5), it follows that ${\displaystyle B\left(h_2\right)<0⟺\mathsf{\Gamma }<\frac{\omega {\left(\mu r-\beta \eta \right)}^2}{\beta k\left(\beta \delta -\gamma r\right)},}$ so that, by (7), one has that $B\left(h_2\right)<0$ if and only if $\mathsf{\Gamma }<{\mathsf{\Gamma }}_b$. Theses (i) and (ii) then easily follow. □

Therefore, if the parameter $\beta$, related to the removal rate of the pathogen by the PMNs, is within the range described by inequalities (8), then it is the parameter $\mathsf{\Gamma }$, (namely the influx of Vitamin D), that can make the difference between the forward and the backward scenario. In fact, if $\mathsf{\Gamma }$ is above the threshold ${\mathsf{\Gamma }}_b$, then a classical forward scenario is obtained. On the contrary, inadequate levels of Vitamin D supplementation lead to a backward scenario.

We stress that Theorem 3 is in perfect agreement with the existence results provided in Theorem 1.

4. Discussions

To critically discuss the implications of the above results, we choose the qualitative set of parameter values summarized in Table 1 and consider the following three cases according to the chosen level of Vitamin D influx.

Case I-Very low Vitamin D influx: ${\displaystyle \mathsf{\Gamma }=0.001<\frac{1}{4}{\mathsf{\Gamma }}_b}$. In this case, a backward scenario is expected with transcritical threshold $h_2=0.0045$ and saddle-node bifurcation threshold $h_1^{\ast }=-0.00178$. The occurrence of such a phenomenology is depicted in the bifurcation diagram at the Top-Left panel of Figure 1. When the influx of Vitamin D is too low, then the saddle node bifurcation value $h_1^{\ast }$ is negative. This means that we always have bistability in the range $[0,h_2]$ whereas, for $h>h_2$, the only attractor is the pathogen-persistent equilibrium. This is surely the worst situation from the point of view of disease recovery. In fact, if $0<h<h_2$ (the Vitamin D degradation is low), a pathogen-free situation can be obtained only in correspondence of a sufficiently low viral load, namely starting in the basin of attraction of the pathogen-free equilibrium. Differently, for higher degradation values ($h>h_2$), system trajectories will tend towards the pathogen-persistent equilibrium, characterized by a high level of pathogen $x_e$. As Figure 1 clearly shows, in this case, decreasing the value of the parameter h below the transcritical threshold $h_2$ will have only the effect of slightly reducing the level of pathogen and a pathogen-free situation cannot be recovered in this case. Therefore if the level of Vitamin D supplementation is too low ($\mathsf{\Gamma }<\frac{1}{4}{\mathsf{\Gamma }}_b$), no matter how much you can reduce its degradation, a complete healthy situation cannot be recovered. This marked characteristic of irreversibility, specific only of this Case I, allows us to link it to situations in which Long-COVID can occur.

Case II-Low Vitamin D influx: ${\displaystyle \frac{1}{4}{\mathsf{\Gamma }}_b<\mathsf{\Gamma }=0.002<{\mathsf{\Gamma }}_b}$. Also in this case, a backward phenomenology is expected with transcritical threshold $h_2=0.009$ and saddle-node bifurcation threshold $h_1^{\ast }=0.00634$. Such a phenomenology is captured by the bifurcation diagram in the Top-Right panel of Figure 1. In this case, the saddle-node bifurcation threshold $h_1^{\ast }$ is positive. A relevant consequence is therefore that, although starting from a fairly low level of pathogen, raising the value of the parameter h above the threshold $h_2$ makes the level of pathogen to suddenly jump to much higher values since system trajectories move towards the pathogen-persistent equilibrium. In this case, a pathogen-free situation cannot be restored by simply reducing the value of h below $h_2$ but it is instead necessary to further lower the level of the degradation parameter h below the saddle-node bifurcation threshold $h_1^{\ast }$. This characteristic of reversibility sub conditione, specific of this Case II, allows us to associate it to situations in which the course of the disease may exhibit sudden complications or a non-gradual recovery of a health status.

Case III-High Vitamin D influx: $\mathsf{\Gamma }=0.1>{\mathsf{\Gamma }}_b$. In this case, a forward scenario is expected, with transcritical threshold $h_2=0.45$. The related phenomenology is depicted in the bifurcation diagram in the Bottom-Center panel of Figure 1. If the level of Vitamin D influx is sufficiently high ($\mathsf{\Gamma }>{\mathsf{\Gamma }}_b$) and its degradation is sufficiently low ($h<h_2$), any initial level of pathogen $x\left(0\right)>0$ will however result in a pathogen-free situation. If the degradation parameter h is incremented to a value $h>h_2$, system trajectories will tend towards the pathogen persistent equilibrium. No jump is observed in this case but a smooth progression towards disease. In this case, a pathogen-free situation can be restored simply by reducing the degradation parameter h below the transcritical bifurcation threshold $h_2$. This marked characteristic of reversibility, specific only of this Case III, allows us to link it to situations in which the disease progresses slowly and there is a gradual recovery towards health.

We observe that in all the three cases (I, II, III), the disease can evolve towards a situation of high severity. For example, a very high level of Vitamin D supplementation (i.e., a very high $\mathsf{\Gamma }$ value as in Case III) can ensure a gradual increase in the disease as well as a smooth recovery from its serious forms, protecting against possible hysteresis effects (as in case II) or from irreversible situations such as Long-COVID (as in Case I) but it is unable to avoid disease progression towards high levels of pathogen. In fact, it is the greater or lesser degradation of Vitamin D in the body—linked to the h parameter—that affects the possibility to develop the disease in a more or less severe form. As a consequence, both these parameters ($\mathsf{\Gamma }$ and h) can hence be considered as possible control parameters for the evaluation of disease recovery and severity.

In this sense, a further interesting point is to identify those factors that can mostly influence the occurrence of the different scenarios. At this regard, we observe that, fixed the other parameters, the parameter $\omega$, namely the degradation of Vitamin D uniquely due to the pathogen, and the parameter k, namely the activation rate of PMNs by Vitamin D, strongly influence the threshold value ${\mathsf{\Gamma }}_b$ as defined in (7). In particular the higher is $\omega$ and/or the smaller is k, the larger is the threshold value ${\mathsf{\Gamma }}_b$. Therefore, high rates of Vitamin D degradation by the pathogen and/or low rates of PMNs activation by Vitamin D favor the occurrence of backward-type scenarios, hence displaying features either of irreversibility (Case I) or of reversibility sub conditione (Case II).

Furthermore, for sufficiently high levels of Vitamin D influx (high $\mathsf{\Gamma }$ values), i.e., such as to ensure a forward-type phenomenology (Case III), high values of the parameter k increase the value of the transcritical threshold $h_2$, therefore enlarging the interval $[0,h_2]$ where the pathogen-free equilibrium is the only attractor for the system. This helps to make the forward-type phenomenology a even more favourable scenario.

Further Remarks on the Subthreshold Regions within the Backward and the Forward Scenario

As just stressed above, in both Case I and Case II (i.e., low values of Vitamin D influx), because of the occurrence of a backward scenario, a bistability range can be found subthreshold. Namely a range of the degradation parameter h is found such that the pathogen-free equilibrium $P_0$ and the pathogen-persistent equilibrium $P^{\ast }$ are both attractors for the system. In this bistability region, the long-term behavior of the system will strongly depend on the initial conditions. A closer inspection to this bistability region, allows us to ask how variations in Vitamin D influx and in the degradation parameter h can impact the basins of attractions [55] of the stable equilibria $P_0$ and $P^{\ast }$.

Figure 2 refers to Case I where a very low Vitamin D influx is considered, i.e., $\mathsf{\Gamma }=0.001$. It is easy to observe that the basin of attraction of the pathogen-free equilibrium $P_0$ decreases with increasing the value of the degradation parameter h whereas that of the pathogen-persistent equilibrium increases.

Moreover, as shown in Figure 3 and Figure 4, for fixed initial values of Vitamin D (i.e., $D\left(0\right)=0.2;0.6;0.8$) the basin of attraction of the pathogen-free equilibrium is wider in correspondence to higher initial values of Vitamin D (i.e., D(0) = 0.8) and is reduced by increasing the values of the degradation parameter h.

The same qualitative features are observed by looking at the bistability region of Case II, namely when higher values of the parameter $\mathsf{\Gamma }$ are considered. This seems to make the management of the bistability regime “somehow” less critical, since starting with a high initial value of Vitamin D and limiting its degradation can help expand the set of initial conditions for which the system can approach a long-term pathogen-free scenario.

Differently in Case III, because of the occurrence of the forward scenario, no bistability regime can be found subthreshold where the pathogen-free equilibrium turns out to be the only attractor for the system. However, also this latter situation can deserve further investigations since transitory changes in the pathogen level could occur before the long-term pathogen-free equilibrium is achieved. Temporary perturbations might in fact induce a transitory increase in the pathogenicity. In this case, trajectory moves away from the stable pathogen-free equilibrium towards higher pathogen levels before eventually coming back to a pathogen free situation. In this case, the pathogen-free equilibrium is called ’reactive’. Exploring reactivity is extremely important to understand the response of dynamic systems to disturbances. It has been investigated in many different contexts as in population dynamics [56,57], in epidemiology [58] and in reaction diffusion systems where it has been recognized as a necessary condition for spatial pattern formation [59]. Reactivity was first introduced in [60] but it has recently been extended to consider linear transformations of the output variable [61] with applications in the ecological field [62,63].

Here we will refer to the standard concept of reactivity [60] that measures the maximum amplification rate of all possible initial perturbations ${\mathbf{x}}_0$ of magnitude $\parallel {\mathbf{x}}_0\parallel$ to an asymptotically stable equilibrium ${\mathbf{c}}^{\ast }$ of a dynamical system ${\displaystyle \frac{d\mathbf{c}}{dt}}=\mathbf{f}\left(\mathbf{c}\right)$, $\mathbf{c}\left(0\right)={\mathbf{c}}^{\ast }+{\mathbf{x}}_0$. More precisely:

Definition 1. 

Given a nonlinear system ${\displaystyle \frac{d\mathbf{c}}{dt}}=\mathbf{f}\left(\mathbf{c}\right)$, $\mathbf{f}:{\mathbb{R}}^N\to {\mathbb{R}}^N$, $\mathbf{c}\in {\mathbb{R}}^N$ and its linearization ${\displaystyle \frac{d\mathbf{x}}{dt}}=J^{\ast }\mathbf{x}$ around an equilibrium point ${\mathbf{c}}^{\ast }\in {\mathbb{R}}^N$ and an initial value $\mathbf{x}\left(0\right)={\mathbf{x}}_0$ the equilibrium is reactive (resp. nonreactive) if the reactivity, defined as

$$r\left(J^{\ast }\right):=\underset{\parallel {\mathbf{x}}_0\parallel \ne 0}{max}\left[{\left.\left({\displaystyle \frac{1}{\parallel \mathbf{x}\parallel }}{\displaystyle \frac{d\parallel \mathbf{x}\parallel }{dt}}\right)\right|}_{t=0}\right],$$

(13)

is positive (resp. negative).

When $\parallel ·\parallel$ denotes the Euclidean norm, the characterization ${\lambda }_{max}\left(H\left(J^{\ast }\right)\right)$ where $H\left(J^{\ast }\right):=({J^{\ast }}^T+J^{\ast })/2$, represents a simple tool to calculate the above quantity.

As the spectral abscissa $a(·)$ (i.e., the maximum of the real parts of the eigenvalues) of a given matrix A is bounded from above by the largest eigenvalue of its hermitian part, i.e., $a\left(A\right)\le r\left(A\right)$, we have that [64]: (i) if the equilibrium point of a linear or nonlinear dynamical system is unstable, then it is reactive; (ii) if the equilibrium point is stable (globally for a linear system or at least asymptotically for a nonlinear system), then it can be either reactive or nonreactive; (iii) if A is normal, then $a\left(A\right)=r\left(A\right)$. Consequently, all the unstable equilibrium points of the system are reactive. Hence, the reactivity needs to be investigated only for stable equilibrium points with a non-normal Jacobian matrix [64].

We want now to explore the reactivity of the pathogen-free equilibrium $P_0$. According to the choice of parameters in Table 1, we have that $J\left(P_0\right)$ in (11) can be written as

$$J\left(P_0\right)=\left[\begin{array}{ccc}-\frac{9\mathsf{\Gamma }}{40h}+\frac{1}{20}& 0& 0\\ -\frac{\mathsf{\Gamma }}{8h}+\frac{19}{12}& -\frac{6}{5}& \frac{3}{10}\\ -\frac{\mathsf{\Gamma }}{2h}& 0& -h\end{array}\right],$$

with eigenvalues $-6/5$, $-h$, $-\frac{9\mathsf{\Gamma }-2h}{40h}$. It makes sense to study the reactivity of $P_0$ when it is a stable attractor for the system i.e., when $\frac{\mathsf{\Gamma }}{h}>\frac{2}{9}$. It results that

$$H\left(J\left(P_0\right)\right)=\left[\begin{array}{ccc}-\frac{9\mathsf{\Gamma }}{40h}+\frac{1}{20}& -\frac{\mathsf{\Gamma }}{16h}+\frac{19}{24}& -\frac{\mathsf{\Gamma }}{4h}\\ -\frac{\mathsf{\Gamma }}{16h}+\frac{19}{24}& -\frac{6}{5}& \frac{3}{20}\\ -\frac{\mathsf{\Gamma }}{4h}& \frac{3}{20}& -h\end{array}\right],$$

and a numerical study suggests that, for this choice of parameters, the largest eigenvalue is always positive, with this meaning that the equilibrium $P_0$ is always reactive (see Figure 5).

As just observed, the dynamics of recovery after a perturbation is characterized by a sudden transitory increase in the pathogen level, followed by its complete disappearance in the long term behavior. Therefore, to ensure the return of trajectories to a pathogen-free situation, not only a stable pathogen-free equilibrium is required but the absence of further perturbations becomes also essential. This leads to some interesting consequences in terms of control since active monitoring and prompt reactions to counteract sudden rises in the pathogen level could be still required even after a pathogen-free situation has been obtained.

5. Conclusions

The role of Vitamin D in the prevention and treatment of COVID-19 has been recently at the center of an intense scientific debate as testified by the many clinical and experimental studies on the topic. However, to the best of our knowledge, no studies have addressed this intriguing question by the means of a mathematical model. Instead, mathematical modeling can be precious to elucidate the role of Vitamin D in the pathogen control because of its capability to provide possible scenarios linked to the adoption of specific clinical measures. This study wants to offer a contribution in this direction with the help of a qualitative phenomenological mathematical model that can allow a better understanding of the mechanisms upon which Vitamin D can interact with the innate immune system and affect its response against SARS-CoV-2.

In our research, we used analytical investigations to understand how Vitamin D dynamics can interact with the immune response in the context of disease recovery and severity. Specifically, we explored how variations in two factors, Vitamin D influx and Vitamin D degradation, can influence the immune system’s ability to fight the disease. The rate at which Vitamin D is introduced into the system plays a crucial role in enhancing the innate immune response. Higher Vitamin D influx levels can lead to an increased production of antimicrobial peptides and cytokines, which aid in combating infections. Our analyses demonstrate that suitably managing Vitamin D influx can contribute to a more robust immune defense, potentially leading to improved disease recovery.

The rate at which Vitamin D is broken down and metabolized affects its availability to the immune system. Slower Vitamin D degradation allows for a sustained presence of active Vitamin D, which has been shown to positively impact immune function. By manipulating Vitamin D degradation rates, we can potentially enhance the immune response’s effectiveness against the disease.

Our analytical investigations reveal that variations in these parameters can have a direct impact on the behavior of the immune system and, consequently, on the progression and severity of the disease. Specifically, by varying Vitamin D influx, the model provides three different dynamical scenarios characterized by different modalities in the recovery of the state of health from the disease. In particular very high level of Vitamin D supplementation are linked to forward-type phenomenologies that are characterized by a complete reversibility, so that disease progresses slowly and, as well, there is a gradual recovery towards health. Instead, by decreasing the level of Vitamin D supplementation, hysteretic phenomenologies via backward scenario emerge, with features of reversibility sub conditione in the recovery from the disease when the Vitamin D supplementation is quite low. This feature is further changed in irreversibility when the Vitamin D influx becomes very low, being compatible with Long-COVID situations.

This point is worthy of special attention because it turns out to be perfectly in line with a recent clinical study [65] in which a connection between Long-COVID syndrome and low levels of Vitamin D has been established. Until now, a limited number of studies have explored the impact of Vitamin D on the likelihood of Long-COVID development. However, these investigations have frequently yielded contradictory findings. Such discrepancies are largely attributed to the inclusion of diverse patient groups, varying in terms of acute disease severity levels, distinct timeframes for follow-up assessments, and disparate collections of symptoms and signs encompassing different aspects of health. In [65] a controlled study with a homogeneous population and a multidisciplinary evaluation was performed, providing evidences that very low levels of Vitamin D can be associated with the presence of Long-COVID syndrome and a greater risk of developing it. The dynamical scenarios obtained with our mechanistic mathematical model are completely in agreement with this findings and suggest how Long-COVID can be the result of hysteretic phenomenologies that could be avoided—or somehow reduced—by suitably increasing the level of Vitamin D influx.

A further point to think about is that, in all the three dynamical scenarios we found (and therefore also in the more favourable scenario with very high Vitamin D influx), the disease can evolve towards a situation of high severity (i.e., high level of pathogenicity) since it is the greater or lesser degradation of Vitamin D in the body (related to the h parameter) that affects the possibility to develop the disease in a more or less severe form. High values of Vitamin D influx (i.e., high $\mathsf{\Gamma }$ values) would eventually contribute to a smooth recovery, reducing the risk of complications and the possible onset of Long-COVID.

In conclusion, moving within the scientific debate that have puzzled physicians about the therapeutic value of Vitamin D for patients with COVID-19, our study confirms—once again—how this hormone can have important effects on the modulation of our body’s immune system and suggests that both Vitamin D influx and Vitamin D degradation should be accurately monitored in the COVID-19 treatment and prevention. In fact, by considering Vitamin D influx and degradation as control parameters, healthcare practitioners and researchers can potentially tailor interventions to optimize the immune response and improve disease outcomes. This insight could hence have implications for the development of therapeutic strategies that leverage Vitamin D dynamics to mitigate the impact of diseases like COVID-19.

It would seem therefore appropriate to monitor the values of circulating Vitamin D in pre and post-COVID patients and to offer supplementation, in case of deficiency. Vitamin D supplementation could therefore represent a possible preventive strategy in reducing disease complications and the burden of COVID-19 sequelae.

Finally, we are surely aware that the model we analyzed barely touches the breathtaking complexity of the biological system under examination since we focused only on a few mechanisms and considered only some essential interactions. Starting from it, efforts are needed in several different directions—for example including adaptive immunity and immunological memories—to gain even more insights on the intriguing relationship between Vitamin D and immunological system. On this line we hope that this first study could inspire other ones, hence becoming a starting point to further explore the multi-pathway and multi-functional role of Vitamin D, also using different mathematical approaches and more refined mathematical models.

Code Availability

The figures depicting bifurcation diagrams and the basins of attraction have been created with MATLAB (version R2023a). The source code will be available upon request.