Generate two negative binomial distributed random variables with predefined correlation

Assume I have a negative binomial distributed variable X1 with NB(mu=MU1,size=s1) and a negative binomial distributed variable X2 with NB(mu=MU2,size=s2). I fitted a negative binomial regression to estimate Mu's and size's from my data

I can use the rnbinom() function in R to generate random draws from this distribution.

X1model<-rnbinom(n=1000,mu=MU1fitted,size=s1fitted)
X2model<-rnbinom(n=1000,mu=MU2fitted,size=s2fitted)

Those draws are now independent. However how can I draw from those distributions, so that they exhibit a predefined correlation r, which is the correlation I observe between my original data X1,X2.

so that:

cor(X1,X2,method="spearman") = r = cor(X1model,X2model,method="spearman")

-or even better draw from those with any arbitrary preset correlation r