Constructing a Lyapunov function for an ODE system that describes epidemic spreading on scale-free networks

Question

I was recently studying an epidemic spreading model, where two competing viruses spread over a scale-free network. $$ \begin{aligned} \frac{dI_{1,k}(t)}{dt} = - I_{1,k}(t) + \psi_1 k (1-I_{1,k} - I_{2,k}) \Theta_1(t)\\ \frac{dI_{2,k}(t)}{dt} = - I_{2,k}(t) + \psi_2 k (1-I_{1,k} - I_{2,k}) \Theta_2(t), \end{aligned} $$ where $\Theta_1(t) = \frac{\sum_{k'}k'P(k')I_{1,k'}}{\langle k \rangle}$ and $\Theta_2(t) = \frac{\sum_{k'}k'P(k')I_{2,k'}}{\langle k \rangle}$.

Here is the interpretation: $k$ represents the degree of a vertex. There are only finite degrees. $P(k)$ is the portion of the vertices that has degree $k$, hence $\sum_{k'}P(k') =1$. $\langle k \rangle \triangleq \sum_{k'} k'P(k')$ is the average degree. $I_{1,k}$ represents the portion of nodes that are infected by virus $1$ among the nodes with degree $k$. $0 \leq I_{1,k},I_{2,k}$ and $I_{1,k}+I_{2,k}\leq 1$.

The term $-I_{1,k}(t)$ in the ODEs gives the recovery speed. The term $\psi_1 k (1-I_{1,k} - I_{2,k}) \Theta_1(t)$ is the transmission speed.

Checking the steady-state by setting $$ \begin{aligned} - I_{1,k}(t) + \psi_1 k (1-I_{1,k} - I_{2,k}) \Theta_1(t) = 0\\ - I_{2,k}(t) + \psi_2 k (1-I_{1,k} - I_{2,k}) \Theta_2(t) = 0, \end{aligned} $$ gives the equilibrium $(I_{1,k}^*,I_{2,k}^*)= (I_{1,k}^*,0)$ when $\psi_1 > \psi_2$, where $I_{i,k}^*$ satisfies the relation $I_{i,k}^* = \frac{\psi_1 k \Theta_1^*}{1 + \psi_1 k \Theta_1^*}$ for all $k$.

The question is how to show this equilibrium is actually globally stable whenever $0< I_{1,k}(0)$. Similation shows global stability. But I have not found a Lypunov function to show theoretical guarantees.

I have tried several forms of Lyapunov candidates $$ \begin{aligned} V(t)= \sum_{k} \left\{ b_1(k) (I_{1,k} - I_{1,k}^*)^2 + b_2(k) (I_{2,k}-0)^2 \right\} + \Theta_1 - \Theta_1^* - \ln \frac{\Theta_1}{\Theta_1^*} + \Theta_2 \\ V(t)= \sum_{k} \left\{ b_1(k) (I_{1,k}- I_{1,k}^* - \ln \frac{I_{1,k}}{I_{1,k}^*}) + b_2(k) I_{2,k} \right\} + \Theta_1 - \Theta_1^* - \ln \frac{\Theta_1}{\Theta_1^*} + \Theta_2 \end{aligned}. $$

I was trying to find $b_1(k),b_2(k)$ which are constant functions of $k$. I tried $b_1(k) = \frac{kP(k)\psi_1}{\langle k \rangle}$, $b_1(k) = \frac{kP(k)\psi_1}{\langle k \rangle}$. Just couldn't find my way to show that $\dot{V} \leq 0$.

Do anybody have any idea about how to find a Lyapunov function for this kind of ODE systems? Or any form of Lyapunov candidates should I look after?

Answer 1

Figured it out after some tries several months ago. I post the answer here in case someone may be interested.

Consider the following Lyapunov function $V(t)$ for $t\geq 0$, $$ \begin{aligned} V(t) =& \frac{1}{2}\sum_{k\in\mathcal{K}} \left[ \tilde{b}_k (S_k(t) - S_k^*)^2\right] + \Theta_1(t) - \Theta_1^* \\ &- \Theta_1^* \ln \frac{\Theta_1(t)}{\Theta^*_1}+\Theta_2(t), \end{aligned} $$ where the coefficients $\tilde{b}_k >0$ are given by $\tilde{b}_k = \frac{kP(k)}{\langle k \rangle S^*_k}$.

Note that $$ \begin{aligned} \dot{\Theta}_1(t) = \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\frac{dI_{1,k}(t)}{dt}=\Theta_1(t)\left[ \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\psi_1 k S_k(t) - 1\right] \end{aligned} $$ and at steady state, $$ \begin{aligned} I_{1,k}^* &= \psi_1 k S_k^* \Theta_1^*,\\ S_k^* & =1- I_{1,k}^*,\\ 1 &= \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\psi_1 k S_k^*. \end{aligned} $$

Then, $$ \begin{aligned} \frac{dV(t)}{dt} &= \sum_{k\in\mathcal{K}}\left[ \tilde{b}_k (S_k - S^*_k)\frac{dS_k}{dt}\right]+ \frac{\Theta_1 - \Theta_1^*}{\Theta_1}\frac{d\Theta_1}{dt} + \frac{d\Theta_2}{dt}\\ &= \sum_{k\in\mathcal{K}} \tilde{b}_k (S_k-S^*_k) \left[ I_{1,k} + I_{2,k} -\psi_2 k S_k\Theta_2 -\psi_1 k S_k\Theta_1\right] +(\Theta_1 - \Theta_1^*)\left[ \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\psi_1 k S_k - 1 \right]\\ &\ \ \ + \Theta_2 \left[ \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\psi_2 k S_k - 1\right]\\ &= \sum_{k\in\mathcal{K}} \tilde{b}_k (S_k - S_k^*)\big[(I_{1,k}-I_{1,k}^*) + I_{2,k}-\psi_2k(S_k\Theta_2 - S_k^* \Theta_2^*) - \psi_1 k(S_k\Theta_1 - S_k^* \Theta_1^*) \big]\\ &\ \ \ + (\Theta_1-\Theta_1^*) \left[ \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\psi_1 k (S_k-S_k^*) \right] + \Theta_2 \left[ \frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k) k (\psi_2 S_k- \psi_1 S_k^*) \right]\\ &= \sum_{k\in\mathcal{K}} \tilde{b}_k \big[(S_k -S_k^*)(I_{1,k} -I_{1,k}^*) +(S_k -S_k^*)I_{2,k}\\ &\ \ \ \ \ \ \ \ \ \ -\psi_2k \Theta_2(S_k -S_k^*)^2 +\psi_2k S_k^* (S_k-S_k^*)(\Theta_2^* -\Theta_2)- \psi_1 k\Theta_1(S_k -S_k^*)^2\\ &\ \ \ \ \ \ \ \ \ \ + \psi_1 k S_k^* (S_k - S_k^*)(\Theta_1^* -\Theta_1) \big]\\ &\ \ \ +\frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k)\psi_1 k (S_k-S_k^*)(\Theta_1-\Theta_1^*)+\frac{1}{\langle k \rangle}\sum_{k\in \mathcal{K}} k P(k) k (\psi_2 S_k- \psi_1 S_k^*)\Theta_2\\ =&\sum_{k\in\mathcal{K}} \tilde{b}_k \big[ (S_k -S_k^*)(I_{1,k} - I_{1,k}^*) + (S_k - S_k^*)I_{2,k}-\psi_2 k\Theta_2(S_k-S_k^*)^2 - \psi_1 k \Theta_1(S_k -S_k^*)^2 \big]\\ &\ \ \ - (\psi_1 -\psi_2) \Theta_2\left[ \frac{1}{\langle k\rangle}\sum_{k\in\mathcal{K}}k^2 P(k )s_k^* \right] \\ \leq &\sum_{k\in\mathcal{K}}\tilde{b}_k \left[(S_k - S_k^*)( 1- S_k -I_{1,k} + I_{1,k} - 1 +S_k^*) \right]- (\psi_1 -\psi_2) \Theta_2\left[ \frac{1}{\langle k\rangle}\sum_{k\in\mathcal{K}}k^2 P(k )S_k^* \right] \leq 0, \end{aligned} $$ where the last inequality holds if $\psi_1 \geq \psi_2$.