  In the last decade, Python-powered physics simulations ecosystem has been
growing steadily, allowing greater interoperability, and becoming an important
tool in numerical exploration of physical phenomena, particularly in soft
matter systems. Driven by the need for fast and precise numerical integration
in colloidal dynamics, here we formulate the problem of Brownian Dynamics (BD)
in a mathematically consistent formalism of the It\=o calculus, and develop a
Python package to assist numerical computations. We show that, thanks to the
automatic differentiation packages, the classical truncated Taylor-It\=o
integrators can be implemented without the burden of computing the derivatives
of the coefficient functions beforehand. Furthermore, we show how to circumvent
the difficulties of BD simulations such as calculations of the divergence of
the mobility tensor in the diffusion equation and discontinuous trajectories
encountered when working with dynamics on $S^2$ and $SO(3)$. The resulting
Python package, Pychastic, is capable of performing BD simulations including
hydrodynamic interactions at speeds comparable to dedicated implementations in
lower-level programming languages, but with a much simpler end-user interface.
