The key for understanding the previous MSE decomposition is to realise that $\theta$ represents a point variable with a fixed value and $\hat{\theta}$ represents a vector containing samples from the estimator (random variable).
As clarification, this formula does not apply to the case where $\theta$ represents a series of values and $\hat{\theta}$ contains element-wise estimations of these points.
A simple numpy example that verifies the previous formula might be useful for understanding this:
import numpy as np
# theta is the target value (point)
theta = 10
# theta_hat contains estimations of theta (vector of 1000 samples)
theta_hat = np.random.normal(theta,0.2,1000)
mse = np.mean(np.square(theta-theta_hat))
var = np.var(theta_hat)
bias = np.mean(theta_hat) - theta
print("MSE result:", mse)
print("Decomposed result:", var + bias*bias)
Which in my case returns 0.042668... in both cases.