Numerical experiments suggest that
$$\sum_{k=1}^{n-1} k^2\sigma(k)\sigma(n-k) = \frac{n^2}{8}\sigma_3(n) - \frac{4n^3-n^2}{24}\sigma(n).$$
PS. In fact, it directly follows from the quoted Touchard and Ramanujan identities.