You should not expect to get the same learned weights (and resulting hidden unit patterns) as they achieve, as it is not unique (you can always do some matched linear transformation of the learned weights) RHW refer to this indirectly in the paper:

It is of some interest that the system that the system employed its ability to use intermediate values in solving this problem. It could, of course, have found a solution in which the hidden units took on only the values of zero and one. Often it does just that, but in this instance, and many others, there are solutions that use the intermediate values, and the learning system finds them even though it has a bias toward extreme values.

A few salient points:

You state that you wish to use "mean-squared loss". However, your code does not do this. When an Elementwise[LogisticSigmoid] is the last layer, NetTrain will default to using the CrossEntropyLoss (see ref/LossFunction) If you wish to use a different loss function that this default, set the LossFunction -> MeanSquaredLossLayer[] option in NetTrain.

You may need to increase the MaxTrainingRounds option in NetTrain. By default (if nothing else is specified) it will run a maximum of 10^4 batches, but this may not be enough. Try increasing this to converge the results. For example, using the default loss (CrossEntropyLoss), I found that the loss was still converging beyond this default; increasing MaxTrainingRounds -> 10^5 batches resulted in a well-converged results.