How to reduce the chemical reaction system $A+B \rightleftharpoons^{k_+}_{k_-} C$ with transfer using the quasi steady state assumption Suppose we have the reaction
$$A+B \rightleftharpoons^{k_+}_{k_-} C$$
within some reaction $\Omega\subset\mathbb{R}^3$, we assume this region to be well mixed and we denote the concentration of $A$ as $a(t)$, similarly for $B$ and $C$. We also know that reactants and complexes can transfer through the surface of $\Omega$, $\partial\Omega$ at some rate. Assuming well mixed and using the law of mass action we may now write
$$\frac{\text{d}a}{\text{d}t}=\underbrace{k_-c-k_+ab}_{\text{reactions}}-\underbrace{ka}_{\text{transfer through } \partial\Omega} $$
$$\frac{\text{d}b}{\text{d}t}=k_-c-k_+ab-kb $$
$$\frac{\text{d}c}{\text{d}t}=-(k_-c-k_+ab)-kc $$
Where k is some rate constant, I want to assume that the kinetic processes within $\Omega$ happen very fast and the system quickly reaches chemical equilibrium, how can I formally reduce this system of ODE's to describe this, i.e. create dependence on only the ratio of the reaction rate constants, using the quasi steady state assumption.
EDIT:
Do we say: Firstly, let us consider the system without elimination of reactants.
$$\frac{\text{d}a}{\text{d}t}=k_-c-k_+ab$$
$$\frac{\text{d}b}{\text{d}t}=k_-c-k_+ab $$
$$\frac{\text{d}c}{\text{d}t}=-(k_-c-k_+ab)$$
we then assume chemical equilibrium, giving
$$c(t)=\frac{k_+}{k_-}a(t)b(t)$$
Therefore returning to our original equations we find that $(k_-c-k_+ab)=0$ thus we have
$$\frac{\text{d}a}{\text{d}t}=-ka $$
$$\frac{\text{d}b}{\text{d}t}=-kb $$
with the condition that 
$$c(t)=\frac{k_+}{k_-}a(t)b(t)$$
EDIT 2:
Or do we in fact say, if we know the clearance rate from $\Omega$ we can define
$$a_{tot}=a+c$$
$$b_{tot}=b+c$$
with equations (derived from our original ODEs)
$$\frac{\text{d}a_{tot}}{\text{d}t}=-ka_{tot}$$
$$\frac{\text{d}b_{tot}}{\text{d}t}=-kb_{tot}$$
Then by assuming chemical equilibrium we know that
$$c(t)=\frac{k_+}{k_-}a(t)b(t)$$
$$c(t)=\frac{k_+}{k_-}(a_{tot}-c)(b_{tot}-c)$$
EDIT 3:
What then confuses me is if we are assuming EDIT 2 is correct we still find that $(k_-c-k_+ab)=0$, therefore from our original equations we find that
$$\frac{\text{d}c}{\text{d}t}=-kc $$
However if $c(t)=\frac{k_+}{k_-}a(t)b(t)$ then this implies that
$$\frac{\text{d}c}{\text{d}t}= \frac{k_+}{k_-}\frac{\text{d}}{\text{d}t}(ab)=\frac{k_+}{k_-}(\frac{\text{d}a}{\text{d}t}b+a\frac{\text{d}b}{\text{d}t})$$
Therefore
$$\frac{\text{d}c}{\text{d}t}=\frac{k_+}{k_-}((-ka)b+a(-kb))=-2k\frac{k_+}{k_-}ab=-2kc\neq -kc$$
What have I done wrong here?
 A: If I understand your question correctly, you want to simplify the ODE system using your knowledge of the time scale difference associated to the processes described by the ODE system. In particular, I understand that you assume that the `internal' reaction dynamics $A + B \leftrightharpoons^{k_+}_{k_-} C$ in $\Omega$ occur on a much faster timescale than the diffusion or leaking of $A,B,C$ through the membrane $\delta \Omega$, i.e. $A \overset{k}{\to} \emptyset $. This would mean that the rate associated with the leaking through the membrane, $k$, is much smaller than both reaction rates $k_{\pm}$, i.e. $k \ll k_{\pm}$. For convenience, I'll rescale the time $t$ with $k_+$, introduce the rate ratios $r=\frac{k_-}{k_+}$ and $\varepsilon = \frac{k}{k_+}$, yielding the system
\begin{align}
 \frac{\text{d} a}{\text{d} t} &= r c - a b &- \varepsilon a,\\
 \frac{\text{d} b}{\text{d} t} &= r c - a b &- \varepsilon b,\\
 \frac{\text{d} c}{\text{d} t} &= -(r c - a b) &- \varepsilon c.
\end{align}
Note that since we assume that $k \ll k_{\pm}$, we have that $\varepsilon \ll 1$.
Using the quasi steady state assumption, we set $\varepsilon$ to zero and look for equilbria of the reduced system. These are given by the solutions of the equation
\begin{equation}
r c = a b.
\end{equation}
To study the effect of the slow diffusion through the membrane, the best thing is to introduce a slow time variable $\tau = \epsilon t$. If we write the original system in terms of this slow time variable, we obtain
\begin{align}
\varepsilon \frac{\text{d} a}{\text{d} \tau} &= r c - a b &- \varepsilon a,\\
\varepsilon \frac{\text{d} b}{\text{d} \tau} &= r c - a b &- \varepsilon b,\\
\varepsilon \frac{\text{d} c}{\text{d} \tau} &= -(r c - a b) &- \varepsilon c.
\end{align}
Again, you see that setting $\varepsilon$ to zero reduces the system to one steady state equation, namely the one we already obtained: $r c = a b$. But what's more important: if we assume we are at (or very close to) the set of equilibria given by $r c = a b$, the resulting system (in slow time $\tau$) reduces to
\begin{align}
 \frac{\text{d} a}{\text{d} \tau} &= - a,\\
 \frac{\text{d} b}{\text{d} \tau} &= - b,\\
 \frac{\text{d} c}{\text{d} \tau} &= - c.
\end{align}
The next observation is that the set of equilibria given by $r c = a b$ is stable. You can show this by linearising the system and determining the signs of the eigenvalues of the resulting $3 \times 3$ matrix, but it's also pretty obvious from the underlying model, which wouldn't really make sense if the reaction equilibrium inside $\Omega$ wouldn't be stable. Anyway, this means that when the system has reached an internal reaction equilibrium in 'fast' time $t$, all components start leaking away slowly through the exponential dynamics described by the slow time ($\tau$) system. This would generally bring the state of the system slowly away from the internal reaction equilibrium $r c = a b$. However, if this perturbation grows too large, the fast internal reaction dynamics 'self-correct' the system towards the internal reaction equilibrium again. To summarise: for any initial concentration $a(0),b(0),c(0)$, the system quickly reaches its 'fast' reaction equilibrium, where $r c = a b$, and then the concentration of every component slowly diminishes, while still obeying the relation $r c = a b$ (up to small perturbations).
Addition: The system can be further analysed as follows: observe that
\begin{equation}
 \frac{\text{d} a}{\text{d} t} + \frac{\text{d} c}{\text{d} t} = -\varepsilon (a + c),
\end{equation}
so $a+c$ decays slowly, and can be solved explicitly as
\begin{equation}
 a + c = (a_0 + c_0)\,e^{-\varepsilon t}.
\end{equation}
Likewise, we have $b + c = (b_0 + c_0)\,e^{-\varepsilon t}$. Hence, we can express both $a$ and $b$ in terms of $c$ and obtain a single ODE for $c$:
\begin{equation}
\frac{\text{d} c}{\text{d} t} = -(r+\varepsilon) c + \left((a_0+c_0) e^{-\varepsilon t} - c\right)\left((b_0+c_0) e^{-\varepsilon t} - c\right).
\end{equation}
You can in theory solve this ODE explicitly, but the solution is rather nasty and doesn't give you too much insight. Instead, we can employ the so-called method of multiple scales. 
This works as follows. The system and our previous analysis gives rise to the idea that $a$, $b$ and $c$ can be most efficiently described using two time scales, $t$ and $\varepsilon t = \tau$ (the 'fast' and 'slow' time scales). The idea is now to treat those two time scales as if they were independent. The smaller $\varepsilon$ is, the more this idea makes sense -- the more the time scales are separated. So, we write $c(t)$ as a function of both $t$ and $\tau$:
\begin{equation}
c(t) = C(t,\varepsilon t) = C(t,\tau).
\end{equation}
Writing $c$ in such a way has impact on the time derivative as well. Using the chain rule, we see that
\begin{equation}
\frac{\text{d} c}{\text{d} t} = \frac{\text{d}}{\text{d} t} C(t,\tau(t)) = \frac{\partial C}{\partial t} + \frac{\text{d} \tau}{\text{d} t} \frac{\partial C}{\partial \tau} = \frac{\partial C}{\partial t} + \varepsilon \frac{\partial C}{\partial \tau}.
\end{equation}
At first, this seems rather counterproductive: we started with an ODE, and now we've got a PDE on our hands! However, this turns out to make our life easier. We can substitute this in the above, and obtain to leading order (i.e. setting $\varepsilon = 0$):
\begin{equation}
 \frac{\partial C}{\partial t} = -r C + \left((a_0+c_0) e^{-\tau} - C\right)\left((b_0+c_0) e^{-\tau} - C\right).
\end{equation}
This differential equation in $t$ has a standard solution of the form
\begin{equation}
 C(t,\tau) = \frac{1}{2}A + \frac{1}{2}\sqrt{B-A^2} \tan\left[\frac{1}{2}\sqrt{B-A^2} (t+\delta(\tau))\right],
\end{equation}
with
\begin{align}
 A &= r+(a_0+b_0+2c_0)e^{-\tau},\\
 B &= 4 (a_0+c_0)(b_0+c_0)e^{-2\tau}.
\end{align}
So, up to an unknown function $\delta(\tau)$, we've found our first order solution for $c$, and in extension for $a$ and $b$.
So, how to find this $\delta(\tau)$? This can be done by using a more detailed assumption for $c$, namely
\begin{equation}
 c(t) = C_0(t,\tau) + \varepsilon C_1(t,\tau).
\end{equation}
You can do exactly the same as above, but now the terms at order $\varepsilon$ give you a new ODE for $C_1$. You can repeat this process ad infinitum, but the equations to be solved will get messier as you go along.
To conclude: it's much easier to stick to the somewhat heuristic description I've given before. Using the exact identities $a + c = (a_0+c_0)e^{-\tau}$ and $b+c = (b_0 + c_0)e^{-\tau}$, you can track the evolution of $a,b,c$ along the set $r c = a b$ by solving $r c = ((a_0+c_0)e^{-\tau} - c)((b_0+c_0)e^{-\tau} - c)$ for $c$ as a function of $\tau$, where in addition $r c_0 = a_0 b_0$. 
As mentioned in the comments: I highly recommend C Kuehn, Multiple Time Scale Dynamics, Springer, 2015, ISBN 978-3-319-12315-8 for further details on this technique and other somewhat equivalent approaches.
