rmk 1. After some manipulation of the series, if
$f(x):=1 - \frac{1-x}{1+x}\sum_{k\in\mathbb Z} x^{k^2}$
Then
$S:=\frac{1}{2}\int_0^1\frac{f(x)}{x} dx.$
So this confirms Charles Matthews' hint:
$\theta:=\sum_{k\in\mathbb Z} x^{k^2}$ Is the theta series, indeed, linked to the Jacobi theta function. (But I don't know how to proceed then.)
rmk 2. Maybe writing Max Muller's series as a sum over double indices (k,n)
$S:=\sum_{1\le k\le n} \frac{1}{(n^2+2k-1)(n^2+2k)}$
one can find a smart partition of the set of indices that allows an iterate summation. (very few chances of success). For instance changing the order of summation, i.e.
$S:=\sum_{k=1}^\infty \sum_{n=k}^\infty \frac{1}{(n^2+2k-1)(n^2+2k)}$
The inner sum can be treated via partial fraction decomposition and summed as illustrated by David Speyer; I found (with $\Psi:=\Gamma'/\Gamma$)
$\sum_{k=1}^\infty\ \left( \mathrm{Im}\frac{\Psi(k+i\sqrt{2k-1})}{\sqrt{2k-1}}-\mathrm{Im}\frac{\Psi(k+i\sqrt{2k})}{\sqrt{2k}}\right).$$\sum_{k=1}^\infty\ \left( \mathrm{Im}\frac{\Psi(k+i\sqrt{2k-1})}{\sqrt{2k-1}}-\mathrm{Im}\frac{\Psi(k+i\sqrt{2k})}{\sqrt{2k}}\right)=0.2666107..$