I am trying to determine the global stability of the following system
$\dot{x} = -r(x+c)(x+y+z)$
$\dot{y} = b_1(a_1 x + (1-a_1)(z-y))$
$\dot{z} = b_2(a_2 x + (1-a_2)(y-z))$
where $b_1,b_2,r\in(0,+\infty)$ and $a_1,a_2,c\in(0,1)$ are fixed parameters. $x$ is bounded in $(-c,+\infty)$
The system has a unique equilibrium point at the origin. My objective is to determine the region in the parameter space for which the origin is globally stable. If I use the following Lyapunov function
$V = \frac{1}{2r}x^2 + \frac{1}{2(b_1(1-a_1)-b_2(1-a_2))}(y+z)^2$
then the derivative turns out to be
$\dot{V} = -\left( y + \frac{x}{2}\left( x+c - \frac{b_1a_1+b_2a_2}{b_1(1-a_1)-b_2(1-a_2)} \right) \right)^2 +\left( z - \frac{x}{2}\left( x+c - \frac{b_1a_1+b_2a_2}{b_1(1-a_1)-b_2(1-a_2)} \right) \right)^2 - x^2(x+c)$
which is obviously not negative definite. I have tried various Lyapunov functions, all of them giving similar expressions for the derivative. For me, even information on boundedness of the solutions of this system would be valuable at this point. If anyone can point me in the right direction, I would be highly grateful.