Best regression model for data that shows exponential decay, but has and is bounded by zero Fairly new to regression analyses. I have a data set for a clinical metric that is bounded in the range of [0-38] for a number of anonymized patients. 
I want to develop a regression model that simulates this clinical metric, but I'm struggling to find the right model given that visually the model follows an exponential decline (e.g., falls to 0 rapidly but in a non-linear way as shown in the picture - each point may represent multiple patients as the points are just superimposed) yet I cannot use a logarithmic/exponential model given that a large number of patients have clinical metric scores of 0. 
Any and all help would be greatly appreciated. 
Thank you!

 A: Here is an implementation in R.  We'll use Poisson regression, though this may not be appropriate for your data. You'll want to do a deviance goodness of fit test on your model once it is fit.  If you reject the null of this test, you should then try a negative binomial regression, or if there are a lot of zeros, a zero-inflated Poisson.  I can add details as necessary.
First, we'll simulate some data. Shown below is the fake data I made
Since these are counts, we'll hypothesize that the distribution of the data on any given data is Poisson with mean $\lambda(t)$. So imagine that the data above are the data you have.  Fit the Poisson regression to that data and continue.  You don't seem to be interested in inference, so I will skip right to the simulation bit.
With the model I fit, I can now predict the mean of the Poisson distribution for any day (remember, the data for any given day is Poisson distributed).  What I'll do is predict those means, then simulate data from Poisson distributions with those means.  Here is the result

The red line is the mean of the distributions, the blue dots are the simulated data.  The intensity of the blue color is a proxy for how frequent the simulated points are.  For instance, at $t = 0$ points are very frequent at $y=10$ but not so frequent at $y = 20$.
Here is the code to make this example.
library(tidyverse)

#Simulate Data for regression
set.seed(1)
times = c(0,sort(sample(0:365, size =  7)))
B = c(log(10),-0.01)
lin.pred = model.matrix(~times)%*%B
y = map(lin.pred, ~rpois(20,exp(.x)))
d = data_frame(times = times, y = y) %>% 
    unnest()

#Plot data
d %>% 
  ggplot(aes(times,y)) + 
  geom_point()

#Fit the poisson regression
model = glm(y~times, data = d, family = 'poisson') 


#Make some simulated data

predict.at = seq(0,365,30)
means = exp(predict(model, newdata = list(times = predict.at)))
simulated_data = map(means, ~rpois(500,.x))
simd = data_frame(times = predict.at, 
                  y = simulated_data) %>% unnest()

lined = data_frame(times = predict.at, y = means)


#Plot the simulated data
simd %>% 
ggplot(aes(times,y))+
  geom_point(color = 'blue', alpha = 0.15)+
  geom_line(data = lined, color = 'red')

Obviously, this should just be a starting point, not an all out solution.  I can't say anything of substance about your data without having it in my hands.  Is this what you were thinking of?
A: It is evident from your graph that your response variable is count data (i.e., non-negative integers).  When you have this kind of response variable it is usual to employ models designed for this type of data, such as a GLM using a negative-binomial distribution or a quasi-Poisson distribution.  A good first step would be to read some materials on regression for count data (see e.g., here, here and here).  In my experience, when dealing with regression data with a count variable as the response, just fitting a standard negative-binomial GLM will usually most of the way to giving you a good model.  In some cases you need to add zero-inflation or some other variation to the model, but it is rare to encounter situations where count data can't be represented with these standard model forms.
