Statistical model of ligand substitution

Recently, I was told that in case of a particular step of a generic ligand substitution reaction:

$$\ce{M(OH2)_{N - n}L_{n} + L <=> M(OH2)_{N - n - 1}L_{n + 1} + H2O}$$

The probability of the forward reaction and by extension, the equilibrium constant of this step, $$K_n$$ would be proportional to

$$\frac{N - n}{n + 1}$$

by a purely statistical analysis. Now I have thought about this for quite a bit, but I can't understand the mathematical reasoning behind arriving at this expression. I suspect it has something to do with the numbers of the ligand being replaced and the ligand which is replacing the other one. Can anyone explain the process of arriving at this expression using simple (if possible) reasoning?

I think the answer may just come down to a simple counting of available sites. The equilibrium constant for a reaction is equal to the ratio of the forward and reverse reaction rates. For the forward reaction, there are $$N-n$$ sites available at which a ligand can replace an $$\ce{H2O}$$. Conversely, for the reverse reaction, there are $$n+1$$ ligand sites at which a water molecules can replace it. If we assume that in each case the reaction rate with $$m$$ sites available is equal to $$m$$ times the reaction rate with $$1$$ site available, we obtain an equilibrium constant that is proportional to the ratio of sites available for the forward and reverse reactions.

$$K=\frac{k_{f,N-n}}{k_{r,n+1}}\approx\frac{(N-n)k_{f,1}}{(n+1)k_{r,1}}\propto\frac{(N-n)}{(n+1)}$$

I think the easy way out is to invoke $$S_\mathrm m = R \ln \Omega$$. If we assume that for a generic complex $$\ce{MA_{n}B_{N-n}}$$,

$$\Omega = {N \choose n} = \frac{N!}{n!(N-n)!} \quad \left[ = {N \choose N-n} \right]$$

and that for the individual molecules $$\ce{A}$$ and $$\ce{B}$$, $$\Omega = 1$$, then the equilibrium constant $$K$$ for

$$\ce{MA_{n}B_{N-n} + A <=> MA_{n + 1}B_{N-n-1} + B}$$

is given by

\begin{align} K &= \exp\left(\frac{-\Delta_\mathrm r G}{RT}\right) \\ &= \exp\left(\frac{\Delta_\mathrm r S}{R}\right) \\ &= \exp\left(\frac{S_\mathrm{m}(\ce{MA_{n + 1}B_{N-n-1}}) + S_\mathrm{m}(\ce{B}) - S_\mathrm{m}(\ce{MA_{n}B_{N-n}}) - S_\mathrm{m}(\ce{A})}{R}\right) \\ &= \exp[\ln\Omega(\ce{MA_{n + 1}B_{N-n-1}}) + \ln\Omega(\ce{B}) - \ln\Omega(\ce{MA_{n}B_{N-n}}) - \ln\Omega(\ce{A})] \\ &= \exp\left[\ln\left(\frac{\Omega(\ce{MA_{n + 1}B_{N-n-1}})\Omega(\ce{B})}{\Omega(\ce{MA_{n}B_{N-n}})\Omega(\ce{A})}\right)\right] \\ &= \frac{\Omega(\ce{MA_{n + 1}B_{N-n-1}})\Omega(\ce{B})}{\Omega(\ce{MA_{n}B_{N-n}})\Omega(\ce{A})} \\ &= \frac{N!}{(n+1)!(N-n-1)!} \cdot \frac{n!(N-n)!}{N!} \\ &= \frac{N-n}{n+1} \end{align}

The reason for ignoring $$\Delta_\mathrm r H$$ is because we are only interested in statistical effects, i.e. entropy, and we don't care about the actual stability of the complex or the strength of the M–L bonds. However, the exact justification for assuming this form for $$\Omega$$ still eludes me. It makes intuitive sense (that there are $$N!/(n!(N-n)!)$$ ways to arrange $$n$$ different ligands in $$N$$ different coordination sites), but I can't convince myself (and don't want to attempt to convince you) that it's entirely rigorous. In particular, I feel like symmetry should play a role here; maybe it is simply that the effects of any symmetry eventually cancel out.