How to determine the distribution of a parameter fit by nonlinear regression 
The example above shows enzyme kinetics -- enzyme velocity as a function of substrate concentration. The well-established Michaelis-Menten equation is:
$Y=V_{max} \cdot \dfrac{X}{K_m + X}$

*

*$X$ are the concentrations of substrate (set by the experimenter)

*$Y$ are the enzyme activities (measured by the experimenter)

*$V_{max}$ is the maximum enzyme velocity at high substrate concentrations. It is fit by nonlinear regression. It has the same units as Y and must be positive.

*$K_m$ is the Michaelis-Menten constant, which is the substrate concentration that leads to half-maximal velocity. Since it is a concentration, it must be positive. It is fit by nonlinear regression, and has the same units as X.

The left panel shows one simulated data set. $V_{max}$ was set to 84 and $K_m$ was set to 4. Each $Y$ value was computed from the equation above plus a random error (Gaussian, SD=12). I made the SD high to make the variation in $K_m$ pronounced. The curve was fit by nonlinear regression using the equation above to determine the $V_{max}$ and $K_m$. Since the residuals are assumed to be Gaussian (and for this example were simulated that way), the nonlinear regression minimizes the sum of the squared residuals.
The middle panel shows the values of $K_m$ fit by nonlinear regression for 100 such simulations. The asymmetry is clear.
The right panel shows the frequency distribution of $K_m$ determined from 10,000 simulations. The distribution was fit to both a normal distribution (red; fits poorly) and a lognormal distribution (blue; fits well). I think this demonstrates pretty clearly that the distribution of $K_m$ is lognormal, or at least it follows a distribution very similar to the lognormal distribution.
My questions are:

*

*For this example, can algebra and/or calculus prove that the distribution of $K_m$ values is lognormal (or prove it has some other distribution)?

*More generally, what method can be used to derive the distribution of any parameter fit by nonlinear regression?

 A: This answer does not (yet) answer the question but should at least help to clarify what the question really is:

"fit by nonlinear regression" sounds like you are using the following model:
$\mathcal{Y}\sim \mathcal{N}(\mu=\frac{X}{X+K_m}, \sigma^2)$

 (this assumes that there is no error in measuring the substrate concentration X; If this nevertheless a good model is another question)

The corresponding likelihood function given a sample $Y^N$ is:
$p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N) = \prod_{i=1}^Np_{\mathcal{N}}(Y^N|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$,
where $p_\mathcal{N}$ is the density of the normal.
and sounds like you are using maximum-likelihood to estimate $K_m$ (and $\sigma^2$).

 (if this is a good approach is yet another question)

$ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N)$
You then seem to sample $\mathcal{Y^N}$ for some fixed $X^N$, $K_m$ and $\sigma$

 (Where  $X^N$ might be your data while $K_m$ and $\sigma$ might be the estimate you obtained for your data with above ML approach)

and then apply above ML estimator ( let's call it ), thus sampling a random variable $\mathcal{\hat{K_m}}$ whose distribution you are asking about (and which you are plotting). There are legit reasons to desire an explicit form of this distribution; for example, to construct confidence intervals for your estimation of $K_m$.

 However since this distribution is not (symmetric and uni-modal) it's yet another question which is the best way to construct a confidence interval given this distribution

Note, however, that this distribution is NOT the posterior distribution of nor a likelihood function for $K_m$ and thus probably not what you desired when you said "the distribution of a parameter".

 the likelihood function is trivial to obtain (look at logLik for your model in R) while the posterior requires you to choose a prior (the empirical distribution of $K_m$ values in databases might be a good choice)

Anyway, let's see how far we get. Let's start by expressing it as compound distribution using the the distribution of $Y^N$ that we know:
$p_{\mathcal{\hat{K_m}}} (\hat{K_M})=\int_{ \{Y^N|\hat{K_M}=ML_{\hat{K_m}}(X^N,Y^N)\}} p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$
This contains $ML_{\hat{K_m}}(X^N,Y^N)$ for which we might be able to find and algebraic expression for:
$ML_{\hat{K_m}}(X^N,Y^N) =  \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \prod_{i=1}^Np_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$
$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(p_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2))$
$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(\frac{1}{\sqrt{2\pi\sigma^2}}) - \frac{\left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2}{2\sigma^2}$
$ = \operatorname*{argmin}\limits_{K_m} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2$
$ 0 = \left.\frac{\mathrm{d}}{\mathrm{d} K_m} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2\right|_\hat{K_m}$
$ = \sum_{i=1}^N \left.\frac{\mathrm{d}}{\mathrm{d} K_m} \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2\right|_\hat{K_m}$
$ = \sum_{i=1}^N  \frac{X^N_i(\hat{K_m}Y^N_i+X^N_i(Y^N_i-1))}{(\hat{K_m}+X^N_i)^3}$
From where I don't know how to continue.

I'm still in the progress of refining this answer please find below a current draft to decide if it's worth your bounty:
In this answer I assume $V_{max}$ is known to be (without loss of generality) 1. As confirmed in the comments you are using the following model:
$\mathcal{Y}\sim \mathcal{N}(\mu=\frac{X}{X+K_m}, \sigma^2)$
The corresponding likelihood function is
$L(K_m, \sigma) = p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N) = \prod_{i=1}^Np_{\mathcal{N}}(Y^N|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$,
where $p_\mathcal{N}$ is the density of the normal distribution.
Now, you would like to know the distribution of a random variable $\mathcal{\hat{K_m}}$ that is the maximum likelihood estimate,
$ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} p_{\mathcal{Y^N}}(Y^N|K_m, \sigma, X^N)$
$ =  \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \prod_{i=1}^Np_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)$
$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(p_{\mathcal{N}}(Y^N_i|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2))$
$ = \operatorname*{argmax}\limits_{K_m} \operatorname*{max}\limits_{\sigma} \sum_{i=1}^N\log(\frac{1}{\sqrt{2\pi\sigma^2}}) - \frac{\left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2}{2\sigma^2}$
$ = \operatorname*{argmin}\limits_{K_m} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+K_m}\right)^2$,
obtained for draws of draws of size $N$ from $\mathcal{Y}$, $\mathcal{Y^N}$, for any $N$, $X^N$, $\sigma$.
You then sampled $K_m$ for some fixed $K$, $X^N$, $K_m$ and $\sigma$ by first sampling $\mathcal{Y^N}$ accordingly and then applying above ML estimator.
Based on this, you think that $\mathcal{K_m}$ follows a log normal distribution.
It is known that, for any differentiable function $f: \mathbb{R}^N \to \mathbb{R}$ and $\mathcal{Y} = f(\mathcal{X})$,
$p_\mathcal{Y}(y) = \int_x \delta(f(x)-y) p_\mathcal{X}(x)\mathrm{d}x$ , where $\delta$ is the Dirac delta.
And that for any monotonic function $g: \mathbb{R} \to \mathbb{R}$ and $\mathcal{Y} = f(\mathcal{X})$,
$p_\mathcal{Y}(y) = p_\mathcal{X}(g^{-1}(y)) \left|\frac{\mathrm{d}}{\mathrm{d}y} g^{-1}(y) \right|$
We can use this to try to derive a closed form for the density of the distribution of $\mathcal{\hat{K_m}}$:
$p_{\mathcal{\hat{K_m}}}(\hat{K_m})=\int \delta (\hat{K_m}-ML_{\hat{K_m}}(X^N,Y^N)) p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$
$\overset{\tiny{\text{if i'm lucky}}}{=}\int \delta(\frac{\mathrm{d}}{\mathrm{d} \hat{K_m}} \sum_{i=1}^N \left(Y^N_i-\frac{X^N_i}{X^N_i+\hat{K_m}}\right)^2) p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$
$=\int \delta(\sum_{i=1}^N  \frac{X^N_i(\hat{K_m}Y^N_i+X^N_i(Y^N_i-1))}{(\hat{K_m}+X^N_i)^3}) p_{\mathcal{Y^N}}(Y^N) \mathrm{d} Y^N$
But I don't how to find a simpler form for that.
For $N=1$ this is a bit simpler:
$p_{\mathcal{\hat{K_m}}}(\hat{K_m})=p_\mathcal{Y}(g^{-1}(\hat{K_m})) \left|\frac{\mathrm{d}}{\mathrm{d}\hat{K_m}} g^{-1}(\hat{K_m}) \right|  =  p_\mathcal{Y}(\frac{X}{X+\hat{K_m}}) \left|\frac{\mathrm{d}}{\mathrm{d}\hat{K_m}} \frac{X}{X+\hat{K_m}} \right|=  p_\mathcal{Y}(\frac{X}{X+\hat{K_m}}) \left|- \frac{X}{(X+\hat{K_m})^2} \right|= p_{\mathcal{N}}(\frac{X}{X+\hat{K_m}}|\mu=\frac{X^N_i}{X^N_i+K_m}, \sigma^2)  \frac{X}{(X+\hat{K_m})^2} $
Where I used:
$ML_{\hat{K_m}}(X^N,Y^N) = \operatorname*{argmin}\limits_{K_m}\left(y-\frac{x}{x+K_m}\right)^2 \Leftrightarrow  0 =\frac{x(\hat{K_m}y+x(y-1))}{(\hat{K_m}+x)^3} \land (\text{further conditions})$ which solves $\hat{K_m}=x(\frac{1}{y}-1)$.
For $N=2$ the explicit form of $ML_{K_m}$ has quite a few more terms
In any case, this shows that $p_{\mathcal{\hat{K_m}}}(\hat{K_m})$ is not log normal (but might converge to it (before converging to normal)).
