where parameters i , i , i (i = 1, 2, 3) and initial conditions x i , y i , z i (i = 0, −1) are nonnegative real numbers. The proposed work is considerably extended and improve some existing results in the literature. Finally, theoretical results are verified by numerical simulations. KEYWORDS boundedness and persistence, equilibrium point, numerical simulations, rate of convergence, stability MSC CLASSIFICATION 40A05; 39A10 ) e −( 3 −2)+ √ ( 3 −2) 2 +4( 1 + 3 ) 2 , Math Meth Appl Sci. 2019;42:7243-7258. wileyonlinelibrary.com/journal/mma