Abstract
This work has considered a mathematical model describing the spread of COVID-19 in a given population. The model comprised 5 systems of equations that take into account different classes describing the impact of COVID-19 in a given population. The time differential operator was replaced with three different types of nonlocal operators. These operators are defined as the convolution of variable order fractal differential operators with different kernels including power law, exponential decay law, and Mittag-Leffler functions. We presented the well-poseness of the models for different differential operators that were presented in detail. A novel numerical scheme was used to solve numerically the system and numerical simulations were provided.