Chaos, Solitons and Fractals 122 (2019) 89–101 Contents lists available at ScienceDirect Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos Behavioural study of symbiosis dynamics via the Caputo and Atangana–Baleanu fractional derivatives Kolade M. Owolabi a,b a Institute for Groundwater Studies, Faculty of Natural and Agricultural Sciences, University of the Free State, Bloemfontein 9300, South Africa b Department of Mathematical Sciences, Federal University of Technology, Akure, Ondo State PMB 704, Nigeria a r t i c l e i n f o Article history: Received 12 January 2019 Revised 2 March 2019 Accepted 17 March 2019 MSC: 34A34 35A05 35K57 65L05 65M06 93C10 Keywords: Atangana–Baleanu–Caputo derivative Chaotic patterns Fractional reaction-diffusion Stability analysis a b s t r a c t Research findings have shown that evolution equations containing non-integer order derivatives can lead to some useful dynamical systems which can be used to describe important physical scenarios. This paper deals with numerical simulations of multicomponent symbiosis systems, such as the parasitic predator- prey model, the commensalism system, and the mutualism case. In such models, we replace the classical time derivative with either the Caputo fractional derivative or the Atangana-Baleanu fractional derivative in the sense of Caputo. To guide in the correct choice of parameters, we report the models linear stability analysis. Numerical examples and results obtained for different instances of fractional power α are pro- vided for non-spatial models as well as the spatial case in one and two dimensions in other to justify our theoretical findings which include the chaotic phenomena, spatiotemporal and oscillatory patterns, mul- tiple steady states and other spatial pattern processes. This paper further suggest an alternative approach to incessant killing of wildlife animals for pattern generation and decorative purposes. © 2019 Elsevier Ltd. All rights reserved. 1. Introduction The study of symbiosis system is not new in literature. It has generated a lot of attentions over the last few decades. Many re- searchers have devoted their attention to the study of popula- tion dynamics which deal with biological or ecological interac- tion of different species of animals or organisms living together in the same habitat. The study of symbiosis models here involve the predator-prey model which describes the situation in which the existence of the predator-species depend solely on the class referred to as the prey, which serve as food to the predatory class. This type of dynamics has been studied widely and had since gen- erated a lot of research interest in areas of applied sciences and engineering. Another class of symbiosis model is called the com- mensalism, in which the existence of one species does not have negative effect on the other. For instance, a typical case describes the relationship between egrets and cattle. The tick (parasitic ani- mal) found on the skin of cows and other livestock serves as food for the egrets. Indirectly, the cows get relieved whenever egrets feed on the ticks, no harm is done to the host. The third example E-mail addresses: mkowolax@yahoo.com, kmowolabi@futa.edu.ng is the mutualism model between two individuals or organisms liv- ing in a given habitat. Consideration is given to mutual love affairs between two species in this work, we refer readers to [1,13,17] for historic details. The aim of this paper is to numerically explore the dynamic richness of multi-species of some symbiosis models such as the predator-prey, commensalism and mutualism systems in fractional medium. The general multispecies symbiosis model takes the form du i dt = F i (u, t ), i = 1, 2, . . . , n (1.1) where u i (t) denotes the population size of the ith species, n stands for the total number of species, and f i accounts for the local kinet- ics. By involving redistribution of species population as a result of random motion of species, the above system of ordinary differen- tial equations become a reaction-diffusion system ∂ u i ∂ t = D i ∂ 2 u i ∂ x 2 + F i (u 1 , u 2 , . . . , u n ), i = 1, 2, . . . , n. (1.2) If the inequality https://doi.org/10.1016/j.chaos.2019.03.014 0960-0779/© 2019 Elsevier Ltd. All rights reserved.