Active Brownian Motion
Coordinates as a matrix
This:
x_co = [np.random.uniform(0, 2) for i in range(n_particles)]
y_co = [np.random.uniform(0, 2) for i in range(n_particles)]
return x_co, y_co

is a little baffling. You call into Numpy to generate uniformly distributed numbers (good) but then put it in a list comprehension (not good) and separate x and y into separate lists (not good).
Consider instead
return np.random.uniform(low=0, high=2, size=(2, n_particles))

and similarly for set_initial_velocities:
return np.random.uniform(low=-1, high=1, size=(2, n_particles))

All of the rest of your code, such as the calculations for xdata/ydata should similarly vectorize over the new length-2 dimension eliminating distinct code for x and y.