Electronic – Understanding the negative feedback and stablizing of a current limiter

Here's a current limiter circuit that uses the negative feedback to limit the current to a certain point.

^{simulate this circuit – Schematic created using CircuitLab}

I have not studied feedback systems before as I'm not an EE student just do it for fun. But here's what I know about this circuit:

When the current through R2 increases the voltage on the base-emitter would also increase. When the voltage reaches about 600mV Q1 will start to conduct. The collector current in Q1 will pull the gate voltage down due to R1 turning the mosfet off, hence reducing the current through the LED. Probably this will continue until Vbe and Vgs reach a certain set point.

Using the equation for a negative feedback loop, assuming the open loop voltage gain is G, we get for the closed-looped voltage gain

$$GV=\frac{G}{1+\beta G}$$
where the denominator takes account for the feedback fraction. Knowing all of this I want to calculate the Vbe and from that measure the current through the LED.

I know the open loop gain of the CE stage but don't know how to start off. I also don't know how to calculate the feedback fraction. I have not seen a practical example of such a circuit before. I want to go through each feedback loop and calculate everything, and finally get a good estimation for all circuit parameters. So any hint would be appreciated.

Here are some observations(out of breadboard):

For Vcc=5V and LED's FWD equal to 2V, I got the following readings:

Vbe=0.47V
Vgs=3.8V
I(led)=40-50mA

UPDATE 1:

I have written a simple code in matlab that runs through a loop for 5000 times, each time calculating a Vbe based on a new Ic.
I used the fact that for a silicon NPN BJT the collector current is 4mA at Vbe=0.7V. For Is I took 8.11*10^-15A. I got a wide range of results as you see below. To save space I've kept some interesting ones. Here's the code and the output:

function clcm()

I=4E-3;
Vbe=0.7;
Ic=40E-6;
Is=8.11E-15;
Vt=26E-3;

R=zeros(10,3);

for i=1:5000
    Vbe=Vbe-Vt*log(I/Ic);
    I=Ic;
    Ic=Is*exp(Vbe/Vt);
    R(i,1)=i;R(i,2)=Vbe;R(i,3)=Ic;
end

disp(num2str(R))

end

Here's the output with the 2nd column indicating the Vbe and the 3rd showing corresponding Ic:

1     0.58026558  3.9954399e-05
2     0.58023592  3.9908849e-05
3     0.58020626  3.9863352e-05
4      0.5801766  3.9817906e-05
5     0.58014694  3.9772512e-05
6     0.58011729   3.972717e-05
7     0.58008763   3.968188e-05
8     0.58005797  3.9636641e-05
9     0.58002831  3.9591454e-05
10     0.57999865  3.9546318e-05

3630      0.4726373  6.3647499e-07
3631     0.47260765  6.3574938e-07
3632     0.47257799  6.3502461e-07
3633     0.47254833  6.3430065e-07
3634     0.47251867  6.3357753e-07
3635     0.47248901  6.3285523e-07
3636     0.47245936  6.3213375e-07
3637      0.4724297   6.314131e-07
3638     0.47240004  6.3069326e-07
3639     0.47237038  6.2997425e-07

4991     0.43227299   1.347568e-07
4992     0.43224334  1.3460317e-07
4993     0.43221368  1.3444972e-07
4994     0.43218402  1.3429644e-07
4995     0.43215436  1.3414334e-07
4996      0.4321247  1.3399041e-07
4997     0.43209505  1.3383766e-07
4998     0.43206539  1.3368508e-07
4999     0.43203573  1.3353267e-07
5000     0.43200607  1.3338044e-07

The Vbe=0.46V is what I exactly measured. I also measured collector current at about 6.8uA so it doesn't agree with what I got above. There should be some constraint that forces me to stop running the code up to a specific iteration, because it's giving me a broad range of currents and voltages. We know the system in a negative feedback should reach equilibrium somewhere.

UPDATE 2:

Here I have plotted the I-V characteristic of the BJT using the data I found above:

Here's the actual behavior of the circuit simulated in LTspice:

where V(n004) is the base-emitter voltage. As you see it doesn't agree with the above schematic. I'm actually trying to understand how I should modify the code to achieve the predicted result.