Numerical Solution to system of balance law with ode system My question is realted to this one, but different. 
Governing equations
I need to simulate the following system, i.e. vector balance law + ode,
$$
\begin{array}{l}
\dfrac{\partial\mathbf{u}}{\partial t} + \dfrac{\partial\mathbf{F}(\mathbf{u})}{\partial x} = \mathbf{G}(\mathbf{u},a,b), \\
\dfrac{da}{dx} = q_1(\mathbf{u},a,b), \\
\dfrac{db}{dx} = q_2(\mathbf{u},a,b). \\
\end{array}
$$
With open boundary conditions (BC) and under some initial conditions (IC).
Questions
Q1.
Suppose, I would like to apply method of lines (MOL) and Runge-Kutta (RK) method. Should I derive differential-algebraic system (DAE), shouldn't I?
Q1.1 Is it possible to discretize initial system on a staggeded (leapfrog) grid by using MOL and RK tecnhiques?
Q2. How to simulate initial system in general?
 A: By open boundary conditions, you maybe mean outflow boundary conditions, which are easy to implement (see e.g. the book by R.J. LeVeque, 2002).
A1. Yes.
A2. Yes, but there may be numerical issues such as a very restrictive stability condition, or numerical stiffness.
A3. You can actually follow the steps in the linked post by introducing a small parameter $\epsilon$ such that
$$
\begin{array}{l}
\dfrac{\partial\mathbf{u}}{\partial t} + \dfrac{\partial\mathbf{F}(\mathbf{u})}{\partial x} = \mathbf{G}(\mathbf{u},a,b), \\
\epsilon\dfrac{\partial a}{\partial t}+\dfrac{\partial a}{\partial x} = q_1(\mathbf{u},a,b), \\
\epsilon\dfrac{\partial b}{\partial t}+\dfrac{\partial b}{\partial x} = q_2(\mathbf{u},a,b). \\
\end{array}
$$
which is a system of balance laws of the form ${\bf q}_t + {\bf \Phi}({\bf q})_x = {\bf R}({\bf q})$. Then, one can split the left-hand side (propagation part) from the right-hand side (relaxation part). For example, one can choose to link the small parameter to the time-step as $\epsilon = \alpha (\Delta t)^\beta$.
