  Probabilistic ordinary differential equation (ODE) solvers have been
introduced over the past decade as uncertainty-aware numerical integrators.
They typically proceed by assuming a functional prior to the ODE solution,
which is then queried on a grid to form a posterior distribution over the ODE
solution. As the queries span the integration interval, the approximate
posterior solution then converges to the true deterministic one. Gaussian ODE
filters, in particular, have enjoyed a lot of attention due to their
computational efficiency, the simplicity of their implementation, as well as
their provable fast convergence rates. In this article, we extend the
methodology to stochastic differential equations (SDEs) and propose a
probabilistic simulator for SDEs. Our approach involves transforming the SDE
into a sequence of random ODEs using piecewise differentiable approximations of
the Brownian motion. We then apply probabilistic ODE solvers to the individual
ODEs, resulting in a pathwise probabilistic solution to the SDE\@. We establish
worst-case strong $1.5$ local and $1.0$ global convergence orders for a
specific instance of our method. We further show how we can marginalise the
Brownian approximations, by incorporating its coefficients as part of the prior
ODE model, allowing for computing exact transition densities under our model.
Finally, we numerically validate the theoretical findings, showcasing
reasonable weak convergence properties in the marginalised version.
