Constructing a Lyapunov function for an ODE system that describes epidemic spreading on scale-free networks 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?
 A: 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$.
