These methods are not really implemented in pytorch so there is no builtin way.
If you have a blackbox optimization method that give you d(loss) / d(out1): how you loss varies wrt the input of Model2, then you can simply do:
grad = BBOpt(out1, out2, loss, ...)
optimizer.zero_grad()
out1.backward(grad)
optimizer.step()

Basically, specify directly what are the gradients that correspond to out1.