My initial answer appears in retrospect to be incorrect and was kindly deleted. This answer resolves the issue, but is considerably slower. The issue arises when NDSolve`EquationProcessing transforms the system of equations into explicit form. Indeed, the problem occurs even, recognizing that once we solve mainEqns for the mainVars, applying that solution to sEqns and zEqns expresses them in explicit form.

The following approaches, based on the knowledge that the last [DoubledGamma] variable and the last mainEqn are special, failed because they took too long (hours in some cases) even for this small example:

1. Eliminate the last [DoubledGamma] from mainEqns using Eliminate. Eliminate takes a very long time.
2. Construct explicit definitions sequentially by a. Solving all but the last mainEqn for all but the last last \[DoubledGamma] variable and then b. Transforming the last mainEqn with that solution and solving the resulting equation for the \[DoubledGamma] variable. That produced explicit expressions in reasonable time, but NDSolve`EquationProcessing takes a very long time to process the definitions even with "SystemSimplification" ->False.
3. Solve with Assumptions, e.g., that z[5,10][t]==0 and z[5,10][t]<s[5,9][t]. With those Assumptions, Solve takes a very long time.

What finally worked is to construct explicit solutions using LinearSolve. That takes, surprisingly, about twice as long as Solve to construct the explicit expressions, but produces the correct ones. To ensure robustness, I added ZeroTest ->(N[#/.(initvals/.Equal->Rule), 60 (*The selected WorkingPrecision*)]==0&) to avoid division by an expression that is 0 under the initial values. Adding that test had little impact on the timing.

My initial answer appears in retrospect to be incorrect and was kindly deleted. This answer resolves the issue, but is considerably slower. The issue arises when NDSolve`EquationProcessing transforms the system of equations into explicit form. Indeed, the problem occurs even, recognizing that once we solve mainEqns for the mainVars, applying that solution to sEqns and zEqns expresses them in explicit form.

The following approaches, based on the knowledge that the last [DoubledGamma] variable and the last mainEqn are special, failed because they took too long (hours in some cases) even for this small example:

1. Eliminate the last [DoubledGamma] from mainEqns using Eliminate. Eliminate takes a very long time.
2. Construct explicit definitions sequentially by a. Solving all but the last mainEqn for all but the last last \[DoubledGamma] variable and then b. Transforming the last mainEqn with that solution and solving the resulting equation for the \[DoubledGamma] variable. That produced explicit expressions in reasonable time, but NDSolve`EquationProcessing takes a very long time to process the definitions even with "SystemSimplification" ->False.
3. Solve with Assumptions, e.g., that z[5,10][t]==0 and z[5,10][t]<s[5,9][t]. With those Assumptions, Solve takes a very long time.

What finally worked is to construct explicit solutions using LinearSolve. That takes, surprisingly, about twice as long as Solve to construct the explicit expressions, but produces the correct ones. To ensure robustness, I added ZeroTest ->(N[#/.(initvals/.Equal->Rule/.exp__[1]:>exp[t]), 60 (*The selected WorkingPrecision*)]==0&) to avoid division by an expression that is 0 under the initial values. Adding that test had little impact on the timing.