Firstly, if a negative binomial model is appropriate, then a model for the second time period with the count the preceding period (ideally log-transformed due to the link function - for 0 counts you may have to e.g. add 0.5) as a model covariate should be fine. This may get quite close to capturing the relationship between preceding and new counts. This model can also be applied when you only look at people with a minimum number of events in the first period. For implementation you can e.g. use proc genmod in SAS or glm.nb from the MASS package in R (just remember that you probably want to correct the standard errors that tend to be too small, confidence intervals and p-values you get from glm.nb). Another option - which is problematic, if you select patients based on the count in the first period, but may be preferable if you did not - is to exploit that a negative binomial model is a Poisson model with a random patient effect. So, you can use a random effects Poisson model with a subject-specific random effect capturing that outcomes for the same patient are correlated across time periods. If you want to stick as closely as possible to a negative binomial model, then a gamma-distributed random effect on the rate would be what you want (e.g. using proc countreg in SAS with the errorcomp=random option in the model statement and the groupid= option in the proc countreg statement). However, you can achieve something very similar with a normally distributed random effect on the log-rate (there is much more software that can do this without much effort, e.g. proc glimmix in SAS or the glmer function in the lme4 package in R).
When evaluating whether a negative binomial model is suitable, remember that standard residuals are not what you want to look at, but look up the recommendations for model residuals for negative binomial models. It also useful to look at the estimated expected proportion with each count from a negative binomial model to to compare it to the data - especially when you are considering whether you need a zero-inflated model.
Secondly, if there is an excess of zero counts and you need a zero-inflated model. Then the same options as above can work. Fitting a zero-inflated model with a covariate is pretty straightforward and most software can do this. A zero-inflated random effects model is a bit trickier, for example you have to decide how the regression coefficients affect the probability of a zero. There seem to be a number of R packages that offer such a model, if you google for “R zero-inflated random effects”, but I do not have any experience with any of them. In SAS I would probably implement this in proc nlmixed by using a normally distributed random effect on the log-event rate, writing out the intended log-likelihood for a zero-inflated Poisson model and using the model ~ general(ll); statement (where ll is the log-likelihood you have defined).
