The study of neural activity synchronization in the brain is a central focus of modern neurobiology and neurodynamics. In a healthy brain, cognitive functions necessitate the accurate integration of neural activity at specific spatiotemporal scales. Some neurological and psychiatric disorders exhibit observable variations in patterns of synchronous population activity. One of the most significant and compelling synchronous population activity patterns is the bursting activity that participates in various informational and physiological processes, including epilepsy. Several mathematical models have been proposed to explain the mechanisms of formation of the bursting activity. The most intriguing and scientifically sound models enable incorporation of astrocytic modulation of neural activity [1]. In this paper, we present a novel phenomenological model for accurately replicating the bursting neural population activity. Our model builds upon the Tsodyks-Markram [2] model and incorporates vital aspects of neuron-glial interaction through a tripartite synapse. The presented model is a simplified version of the previously proposed model [3]. During astrocyte activity (which lasts seconds), we can set the probability of neurotransmitter release, u, which is determined in fractions of a second. This enables the proposed model to be written in the following manner: r=–E+α ln(1+exp((Ju(Y)xE+I0)/α)), =(1–x)/τD–u(Y)xE, =–y/τY+βσy(X). Here, E(t) represents the average neuronal activity of the excitatory population, while x(t) denotes the quantity of available neurotransmitter. Additionally, y(t) describes the concentration of gliotransmitters, which are released as a result of biochemical reactions during neuron-astrocyte interactions. The function describes the alteration in the likelihood of neurotransmitter release when gliotransmitters are present u(y)=u0+(∆u0)/(1+exp(–50(y–ythr)). Here u0 represents the probability of neurotransmitter release without astrocytic influence; ∆u0 indicates the change in the probability of neurotransmitter release caused by the interaction of the gliotransmitter with the presynaptic terminal, and ythr represents the threshold value that determines the change in the probability of neurotransmitter release due to the effect of a gliotransmitter. The function describes the influence of neurotransmitter on the concentration of the gliotransmitter σy(x)=1/(1+exp(–20(x–xthr)), where xthr is the astrocyte activation threshold. Note that the proposed model does not account for synaptic depression’s mechanism. Therefore, various dynamical regimes’ formation in the model is solely controlled by astrocytic dynamics. In the study, control parameters I0 and u0 were selected while the remaining parameters were held constant. Specifically, τ=0.013, τD=0.08, α=1.58, J=3.07. Changes in the concentration of neurotransmitters and gliotransmitters were characterized by parameters ∆u0=0.305, τν=3.3, β=0.3, xxthr=0.75, ythr=0.4. The research was conducted using computer modeling and numerical methods in nonlinear dynamics. The proposed model exhibits a diverse range of population activity patterns, such as spiking and bursting regular and chaotic activity. Chaotic dynamics zones are identified in the control parameter plane. An account is given of the mechanisms behind the emergence of these activity patterns through bifurcation analysis. It is shown that the occurrence of chaotic activity in the system can be associated both with the cascade of period doubling bifurcations and with the further development of chaos, as a result of which a homoclinic attractor appears in the system according to the Shilnikov scenario. It was also shown that multistability is observed in the system in a certain range of parameters. Note that the emergence of multistability and bursting activity in the model is independent of the intricacy of neuron and glial cell dynamics. These dynamics are instead influenced by the feedback loop between the presynaptic neuron and the glial cell. The effects of neuron-like dynamics and neural-glial interaction demonstrated are generally applicable, without specifying particular characteristics of neural-glial interaction, network architecture, or single-cell neuron dynamics. It should be noted that our proposed model solely concerns the potentiation of synaptic transmission by astrocytes. In summary, the proposed phenomenological model for population activity can replicate various patterns of neuron activity found in a wide range of dynamic memory and information processing studies. This research can potentially lead to the development of new, effective treatment methods for neurological diseases that involve neuron-glial interaction. Another potential application for these results is in the development of an optimized live chip that has preset functions. This requires a deeper understanding of rhythmogenesis in neural networks and brain function.