Q:

Python updating numpy array: type error numpy.int32 is not iterable

I have written a function and while loop which is designed to simulate the following system of chemical reactions over time: 
CHEMICAL EQUATIONS IN SYSTEM:
E + S --> ES == 1E + 1S + 0ES + 0P --> 0E + 0S + 1ES + 0P
ES --> E + S == 0E + 0S + 1ES + 0P --> 1E + 1S + 0ES + 0P
ES --> E + P == 0E + 0S + 1ES + 0P --> 1E + 0S + 0ES + 1P
The system is expressed using the following arrays: 
discrete population numbers of Enzyme, Substrate, Enzyme-Substrate complex, Product (at time= 0):
popul_num = np.array([200, 100, 0, 0])

Ratio of reactants: 
LHS = np.array([[1,1,0,0], [0,0,1,0], [0,0,1,0]])

Ratio of products:
RHS = np.matrix([[0,0,1,0], [1,1,0,0], [1,0,0,1]])

Rate of the three reactions in system: 
stoch_rate = np.array([0.0016, 0.0001, 0.1])

state change array:
state_change_array = RHS - LHS

The maximum simulation time: 
tmax = 20
and the start time of the simulation: 
tao = 0.0
I then have the following function and while loop to calculate the how the discrete numbers of molecules in popul_num change over time during the course of the above reactions. 
The popul_num array is then updated with the new values after the reactions have run for a certain period of time: 
def propensity_calc(LHS, popul_num, stoch_rate):
propensity = np.zeros(len(LHS))
for row in range(len(LHS)):
        a = stoch_rate[row]     # type = numpy.float64
        for i in range(len(popul_num)):
            if (popul_num[i] >= LHS[row, i]):       
                binom_rxn = binom(popul_num[i], LHS[row, i])
                a = a*binom_rxn
            else:
                a = 0
                break
        propensity[row] = a     # type = numpy.ndarray
return propensity

propensity = np.zeros(len(LHS))
while tao < tmax:
  propensity = propensity_calc(LHS, popul_num, stoch_rate)
  a0 = (sum(propensity))
  if a0 == 0.0:
      break
  t = np.random.exponential(1/a0)
  rxn_probability = propensity / a0   # propensity = array a0 = number --> Error
  num_rxn = np.arange(rxn_probability.size)
  if tao + t > tmax:
      tao = tmax
      break
  j = stats.rv_discrete(values=(num_rxn, rxn_probability)).rvs()
  print(tao, t)
  tao = tao + t
  popul_num = popul_num + np.squeeze(np.asarray(state_change_array[j]))

I want to use matplotlib to plot a separate line on a line graph for each of the species in the popul_num showing their changing numbers over time. 
I tried appending each of the new values into popul_num and plotting the new array using numpy.append() but didnt have any luck.
I've tried writting a loop to iterate over each element of the array as follows: 
for i in range(4):
plt.plot(list(enumerate(popul_num[i])))  # error numpy.int32 object is not iterable
plt.show()

but I get the following error: 
TypeError: 'numpy.int32' object is not iterable
Is there a way of solving this or a different way to create a plot of popul_num for each species over time?
EDIT: 
At the moment this is the plot that the code is producing 
This is the type of plot that want but i want a separate line for each element of the popul_num array showing the change in that molecular species over time. 
cheers

A:

The popul_num are 4 values that are overwritten in each step. To draw curves, you need to save all these values in a separate array. The following code shows an example:
import matplotlib.pyplot as plt
import numpy as np
from scipy.special import binom
from scipy import stats

popul_num = np.array([200, 100, 0, 0])

LHS = np.array([[1, 1, 0, 0], [0, 0, 1, 0], [0, 0, 1, 0]])
RHS = np.array([[0, 0, 1, 0], [1, 1, 0, 0], [1, 0, 0, 1]])
stoch_rate = np.array([0.0016, 0.0001, 0.1])
state_change_array = RHS - LHS

def propensity_calc(LHS, popul_num, stoch_rate):
    propensity = np.zeros(len(LHS))
    for row in range(len(LHS)):
        a = stoch_rate[row]  # type = numpy.float64
        for i in range(len(popul_num)):
            if (popul_num[i] >= LHS[row, i]):
                binom_rxn = binom(popul_num[i], LHS[row, i])
                a = a * binom_rxn
            else:
                a = 0
                break
        propensity[row] = a  # type = numpy.ndarray
    return propensity

tmax = 20
tao = 0.0

popul_num_all = [popul_num]

propensity = np.zeros(len(LHS))
while tao < tmax:
    propensity = propensity_calc(LHS, popul_num, stoch_rate)
    a0 = (sum(propensity))
    if a0 == 0.0:
        break
    t = np.random.exponential(1 / a0)
    rxn_probability = propensity / a0  # propensity = array a0 = number --> Error
    num_rxn = np.arange(rxn_probability.size)
    if tao + t > tmax:
        tao = tmax
        break
    j = stats.rv_discrete(values=(num_rxn, rxn_probability)).rvs()
    print(tao, t)
    tao = tao + t
    popul_num = popul_num + np.squeeze(np.asarray(state_change_array[j]))
    popul_num_all.append(popul_num)

popul_num_all = np.array(popul_num_all)
for i, (label, color) in enumerate(zip(['Enzyme', 'Substrate', 'Enzyme-Substrate complex', 'Product'],
                                       ['limegreen', 'dodgerblue', 'orange', 'crimson'])):
    plt.plot(popul_num_all[:, i], label=label, color=color)
plt.legend()
plt.tight_layout()
plt.show()

A legend is added with the name of each curve.

