## Appendix A Mathematical Proofs  
### Proof of Theorem 1  
Proof.: First, we have \(\left.\frac{dN}{dt}\right|_{N=0}=0\) from (11). According to (9) and (10), we also have the following:  
\[\left.\frac{dN_{S}}{dt}\right|_{N_{S}=0}=N\max\{C_{w}-C_{g},0\}\geq 0\]  
Since \(N\) is continuous and \(\left.\frac{dN}{dt}\right|_{N=0}=0\), then the solutions \(N\) will remain unchanged once they reach \(N=0\). Hence, for initial conditions \(N(0)>0\), the solutions cannot be negative, and therefore, we conclude that \(N\) is positive invariant. Now, we shall prove positive invariance for \(N_{g}\). For \(N_{g}(0)\), we can have three cases regarding the relative costs of gamergates and non-reproductive workers: (1) \(C_{w}>C_{g}\), (2) \(C_{w}=C_{g}\), or (3) \(C_{w}<C_{g}\). One of these three cases must happen initially, and \(C_{w}=C_{g}\) would be the transition point between cases 1 and 3. Let \(t_{0}\) be the initial time, and \(t_{1}\) be the time when a transition between cases happen. Then, we prove that no matter the case the system starts with, for positive initial conditions the solution \(N_{g}(t)\) is positive.  
1. \(C_{w}>C_{g}\). Let \(\epsilon(t)=C_{w}(t)-C_{g}(t)\), then \[\frac{dN_{g}}{dt}=N_{g}\bigg{[}\bigg{(}\frac{N}{N_{g}}-1\bigg{)}\epsilon-d_{g}N \bigg{]}\] Solving the above differential equation gives \[N_{g}(t)= \exp\bigg{(}\int_{t_{0}}^{t}(-\epsilon(s)-dgN(s))ds\bigg{)}*\] (A1) \[\bigg{[}N_{g}(t_{0})+\int_{t_{0}}^{t}exp\bigg{(}\int_{t_{0}}^{r}( -\epsilon(s)-dgN(s))ds\bigg{)}\epsilon(r)N(r)dr\bigg{]}\] \[\geq 0\]
2. \(C_{w}=C_{g}\), then \[\frac{dN_{g}}{dt}=-d_{g}NN_{g}\] Solving the above differential equation gives \[N_{g}(t)=N_{g}(t_{0})exp\bigg{(}-d_{g}\int_{t_{0}}^{t}N(s)ds\bigg{)}\geq 0\] (A2)
3. \(C_{w}<C_{g}\). Let \(\epsilon(t)=C_{w}(t)-C_{g}(t)\), then \[\frac{dN_{g}}{dt}=N_{g}\big{(}-\epsilon-d_{g}N\big{)}\] Solving the above differential equation gives \[N_{g}(t)=N_{g}(t_{0})\exp\bigg{(}\int_{t_{0}}^{t}(-\epsilon(s)-dgN(s))ds\bigg{)}\geq 0\] (A3)  
The above reasoning proves that no matter which case the system falls in initially, the solution \(N_{g}\) will be positive. Namely,  
\[N_{g}(t)=\left\{\begin{array}{rl}\text{(A1)}&\text{if $t$ such that $C_{w}(t)>C_{g}(t)$}\\ \text{(A2)}&\text{if $t$ such that $C_{w}(t)=C_{g}(t)$}\\ \text{(A3)}&\text{if $t$ such that $C_{w}(t)<C_{g}(t)$}\end{array}\right.\] (A4)Note that at \(C_{w}=C_{g}\), the system has only one equilibrium point \((N^{*},N_{g}^{*})=(0,0)\) which is unstable. Thus, the system will not stay in the case \(C_{w}=C_{g}\), but it will make a transition to either cases 1 or 3. Whichever case the system transitions to, at any time \(t>0\), the solutions will be positive as proved before, and therefore we could conclude that our proposed system (9) and (10) is positively invariant in \(\mathbb{R}^{+2}\).  
Note that \(N\geq N_{g}\), then we have:  
\[\frac{dN}{dt} = eN_{g}\frac{\beta_{n}N^{h_{w}}}{\beta_{g}N_{g}^{g_{g}}+\beta_{n}N ^{h_{w}}}-d_{g}N_{g}N-d_{w}(N-N_{g})N\] \[\leq eN_{g}-d_{g}NN_{g}\] \[\leq N(e-d_{g}N)=Ne\bigg{(}1-\frac{N}{e/dg}\bigg{)}\]  
This implies that  
\[\limsup_{t\to\infty}N(t)\leq\frac{e}{d_{g}}.\]  
Therefore, the colony size, \(N\), is bounded when time is large enough. As a consequence, we could conclude that \(N_{g}\) and \(N_{w}\) are also bounded since \(N=N_{g}+N_{w}\).