Nonlinear Autoregressive model parameter estimation from time series

I'm working on a nonlinear multivariate autoregressive model of order 1 (markovian). It is a discrete-time dynamical system which models exchange of mass between compartments in a compartmental model where each compartment corresponds a given population type characterized, for every time instant, by a fraction of mass $x_u(t)$. The network dynamics are expressed by quantities of mass which have chances of interact and being converted to another type, thus leading to coupled interactions which determine flow of mass between types.

The mass for type $u$ at time $t + 1$ is given by its mass at $t$ plus the incoming flows, coming from a subset of neighboring types, minus outgoing flow to other types.

$$x_u(t+1) = x_u(t) + \sum_{v : u \in N_{\mathcal C}(v)} f_{vu}(t) -\sum_{v \in \mathcal N_{\mathcal C}(u)} f_{uv}(t)$$

Where each flow has the form:

$$f_{uv}(t) = x_u (t) \cdot \sum_{z \in S_u} \frac{x_z(t)}{(\alpha_{u} - 1 )\cdot x_u(t) + B } \cdot \frac{p_{uz}}{m}$$

My objective is doing estimation of parameters $\alpha_u$ (1 for each compartment) and $p_{uv}$ (one for each couple of compartments) from a time series of observations concerning the masses $x_u$. I obtained good results using least squares nonlinear parameters estimation for low dimensional systems, while it quickly becomes impractical for systems with moderate to high number of compartments. Could someone suggest a technique/direction to learn the parameters of this kind of models?

EDIT: the parameters $p_{uv}$ are constrained to be probabilities (thus ranging in [0,1]) and $\alpha_u \geq 1$

Let's denote $\hat{x}_u(t)$ the mass predicted by the model given the parameter estimates $\hat{\alpha}_u$ and $\hat{p}_{uv}$. Then, you can define the following function:
$$S(\alpha,p_{uv};x_u) = \sum_{u,t} \left( x_u(t) - \hat{x}_u(t) \right)^2 \,.$$
This function returns a scalar that is the sum over all $u$ and $t$ of the squared differences between the observed values and the predicted values. Thus, it is a measure of the prediction errors given a set of parameter estimates. (The square is taken to avoid that negative errors are offset positive errors.)