# Update gradients with respect to weights and biases
dWz += np.dot(dz, combined.T)
dWr += np.dot(dr, combined.T)
dWh += np.dot(dh_, combined_r.T)

dbz += dz.sum(axis=1, keepdims=True)
dbr += dr.sum(axis=1, keepdims=True)
dbh += dh_.sum(axis=1, keepdims=True)