Is this multivariate normal? 2 time series linked by a common process

Summary: Consider a scenario where you observe the inputs ($X$) to and outputs ($Y$) from a process ($B$). If I have a model describing how $X$ evolves over time, and a similar model for $Y$, how do I combine the two?

Example + Details To make this concrete, consider the coupled $B$iological process of photosynthesis + respiration ($B = B_p + B_r$), where, to a rough approximation, for every atom of Carbon ($X$) taken up, 1 atom of Oxygen ($Y$) is released.

There are other factors ($F_x$, $F_y$) that also play minor roles in the abundance of $X$ or $Y$, but these "factors" are deterministic processes based on $X$ or $Y$ at the previous time step and 2 known control vectors ($K$, $Z$), and the parameters of interest are associated with the process $B$. $B$ is a function of 2 control variables (in a Tx2 matrix, $U$), and 2 parameters (say, $\beta_1$ and $\beta_2$, in a 2x1 matrix, $\beta$).

The model for $X$ might look something like $$X_t = X_{t-1} + B_{t} + F_{x,t}$$ $$B_{t} = U_t\beta$$ $$F_{x,t}=K_{x,t}(Z_{t-1}-X_{t-1})$$ Similarly, $Y$ is driven by $B$, but in the opposite direction $$Y_t = Y_{t-1} + -B_{t} + F_{y,t}$$$$B_{t} = U_t\beta$$ $$F_{y,t}=K_{y,t}(Z_{t-1}-Y_{t-1})$$

Current Framework, Objective I already have both models for $X$ and $Y$ working in several frameworks (as a dlm fit with maximum likelihood, in a Kalman filter, and in a Bayesian framework). But I want to put the two analyses (for $X$ and $Y$) together.

JAGS Example In JAGS, the main loop of my current model for $X$ looks something like this, and would look identical for $Y$ (the signs of estimated parameters reversed):

for(i in 2:N){
X[i] ~ dnorm(a[i], tauV) # observations (X) are distributed with mean of true values, but with precision tauV (1/tauV is variance of observation error)
Kx[i] ~ dnorm(KxP[i-1, 1], 1/KxP[i-1, 2]) #distributin on Kx, from input KxP containing means and variances of Kx at each time step
Fx[i] <- Kx[i]*(Z[i-1]-X[i-1]) # Z is just an input control vector of know values
aHat[i] <- a[i-1] + U[i,]%*%Beta + Fx[i] # the process
a[i] ~ dnorm(aHat[i], tauW) # true values have a mean of estimated values, but accompanied by process error (process precision is tauW)
}

Motivation It seems that I am not making use of valuable information when I estimate $\beta$ independently from $X$ and $Y$, especially b/c both time series are prone to different sources of observation and process error (as well as some of the same sources, especially for the process error).

Question How do I couple the models? Do I model the observations as multivariate normal?