From microbial communities, human physiology to social and biological/neural networks, complex interdependent systems display multi-scale spatio-temporal patterns that are frequently classified as non-linear, non-Gaussian, non-ergodic, and/or fractal. Distinguishing between the sources of nonlinearity, identifying the nature of fractality (space versus time) and encapsulating the non-Gaussian characteristics into dynamic causal models remains a major challenge for studying complex systems. In this paper, we propose a new mathematical strategy for constructing compact yet accurate models of complex systems dynamics that aim to scrutinize the causal effects and influences by analyzing the statistics of the magnitude increments and the inter-event times of stochastic processes. We derive a framework that enables to incorporate knowledge about the causal dynamics of the magnitude increments and the inter-event times of stochastic processes into a multi-fractional order nonlinear partial differential equation for the probability to find the system in a specific state at one time. Rather than following the current trends in nonlinear system modeling which postulate specific mathematical expressions, this mathematical framework enables us to connect the microscopic dependencies between the magnitude increments and the inter-event times of one stochastic process to other processes and justify the degree of nonlinearity. In addition, the newly presented formalism allows to investigate appropriateness of using multi-fractional order dynamical models for various complex system which was overlooked in the literature. We run extensive experiments on several sets of physiological processes and demonstrate that the derived mathematical models offer superior accuracy over state of the art techniques.