Using Riemann-Stieltjes methods for integrators of bounded $p$-variation we define a pathwise integral driven by a fractional Lévy process (FLP). To explicitly solve general fractional stochastic differential equations (SDEs) we introduce an Ornstein-Uhlenbeck model by a stochastic integral representation, where the driving stochastic process is an FLP. To achieve the convergence of improper integrals, the long-time behavior of FLPs is derived. This is sufficient to define the fractional Lévy-Ornstein-Uhlenbeck process (FLOUP) pathwise as an improper Riemann-Stieltjes integral. We show further that the FLOUP is the unique stationary solution of the corresponding Langevin equation. Furthermore, we calculate the autocovariance function and prove that its increments exhibit long-range dependence. Exploiting the Langevin equation, we consider SDEs driven by FLPs of bounded $p$-variation for $p<2$ and construct solutions using the corresponding FLOUP. Finally, we consider examples of such SDEs, including various state space transforms of the FLOUP and also fractional Lévy-driven Cox-Ingersoll-Ross (CIR) models.