I’m not sure I understand exactly what you want to do, but would the same reparametrisation trick as in the VAE paper and implementations (e.g. pytorch/examples) work with the “usual” procedure?
You would convert standard normal randoms to a variable and then transform them with mu and sigma2. That way, the randoms are fixed w.r.t. the differentiation.