Recently I was studying on the global stability of a model system given as $$\frac{\partial S}{\partial t}=rS\bigg(1-\frac{S+I}{K}\bigg)-\frac{\beta SI}{S+I+c}+D_s\nabla^2 S\\ \frac{\partial I}{\partial t}=\frac{\beta SI}{S+I+c}-aI+D_I\nabla^2I$$ where $p=(x,y)\in \Omega=[0,L]\times [0,L]$ and $t>0$ and $\nabla^2$ being the two-dimensional Laplace operator, with zero-flux boundary conditions $\displaystyle{\frac{\partial S}{\partial n}=\frac{\partial I}{\partial n}=0}$.
Now a positive definite function $V(S,I)$ about the equilibrium point $(S^*,I^*)$ is defined as $$V(S,I)=\int_{S^*}^S\frac{\xi -S^*}{\xi}d\xi+k_1\int_{I^*}^I \frac{\eta-I^*}{\eta}d\eta$$ for some $k_1>0$.
Finally the Lyapunov function is defined by $$E(t)=\iint_\Omega V(S,I)dA$$ on the plane region $A$.
The simplification now goes as follows $$\frac{dE(t)}{dt}=\iint_\Omega \bigg[\frac{\partial V}{\partial S}\frac{\partial S}{\partial t}+\frac{\partial V}{\partial I}\frac{\partial I}{\partial t}\bigg]dA$$
Now afterwards the calculations are becoming too much cumbersome to show that $\displaystyle{\frac{dE(t)}{dt}}$ will be positive or negative to determine the stability. Can someone help me to simplify this derivative as short as possible, by some tricks (if there's any) in multivariable calculus? This will be much helpful. Thanks for your time.