The tests in test_model.py are not working anymore. They should be either updated or removed.

The model should be contained in a class, so that we can easily instantiate independent models in the same script.

When I test the code with one trial, I get very high firing rate of CJ quickly.
At t = 0.002s, the firing rate is already at ~16 spikes/s.
And at t= 0.0055s, the firing rate exceeds 166 spikes/s.
Looking at phi (equation 6), at t = 0.002s, I have phi ~ 57,55.
In the equation, the only variable is I_syn, so I look at it.
All the terms of I_syn are low compared to I_ampa_ext (equation 7).
But I_ampa_ext is determined by a constant term and the white noise, that is negligible.

Would there be some errors in the code units?

The code should have unit tests, so that we can quickly see if a change in the code translate into a change of behavior. This is especially important when we want to do changes that should not effects on the code.

Orientation of axes may be change to fit with those of the article.
Example from the figure 4 :

The cj_b_cells and cj_a_cells function are identical, save for print statements. Only one of them should exist.

Find a way to find the corresponding (quantity of juice A, quantity of juice B) on the abscissa for each firing rate on the ordinate for the 3e column of the figure 4.
Example of abscissa from the figure 4C :

These are the graphs I should obtain with a firing rate that is stable between 3 and 4 spikes/s before the offer and that is under ~25 spikes/s after the offer for CJ cells.

But before the offer I obtain a firing rate between 3 and 4 spikes/s but after the offer I obtain a firing rate that can go until 34 spikes/s.

Concerning CV cells, the firing is stable at ~9,5 spikes/s instead of ~11 spikes/s before the offer and after the offer it can rise until ~16 spikes/s instead of ~18 spikes/s.