diff --git a/BoolODE/__init__.py b/BoolODE/__init__.py
index 70f71a4..f3be3ab 100644
--- a/BoolODE/__init__.py
+++ b/BoolODE/__init__.py
@@ -67,6 +67,7 @@ def __process_jobs(self) -> Dict[int, Dict]:
         Default parameter values are specified here.
         '''
         jobs = {}
+        experimental_settings = {'noise_strength':10}
         for jobid, job in enumerate(self.job_settings.jobs):
             data = {}
             # Create output folder if it doesnt exist
@@ -99,7 +100,15 @@ def __process_jobs(self) -> Dict[int, Dict]:
             data['add_dummy'] = job.get('add_dummy',False)
             data['max_parents'] = job.get('max_parents',1)
             data['modeltype'] = self.global_settings.modeltype
+            data['noise_strength'] = job.get('noise_strength', 10)
 
+            warnlist = [setting for setting in data.keys()\
+                        if (setting in experimental_settings.keys())\
+                        and (data[setting] != experimental_settings[setting])]
+            for ef in warnlist:
+                print(f"___WARNING___\ntIN JOB [{data['name']}] SETTING [{ef}: {data[ef]}]")
+                print(f"={ef}= is currently an experimental feature in BoolODE!")
+                print(f"Results obtained by varying this parameter are\nnot guaranteed to produce optimal results!")
             jobs[jobid] = data
         return(jobs)
 
diff --git a/BoolODE/run_experiment.py b/BoolODE/run_experiment.py
index f40b54a..7a31ef3 100644
--- a/BoolODE/run_experiment.py
+++ b/BoolODE/run_experiment.py
@@ -129,6 +129,7 @@ def Experiment(mg, Model,
     argdict['proteinIndex'] = proteinIndex
     argdict['revvarmapper'] = revvarmapper
     argdict['x_max'] = mg.kineticParameterDefaults['x_max']
+    argdict['noiseStrength'] = settings['noise_strength']
 
     if settings['sample_cells']:
         # pre-define the time points from which a cell will be sampled
@@ -297,6 +298,7 @@ def simulateAndSample(argdict):
     seed = argdict['seed']
     pars = argdict['pars']
     x_max = argdict['x_max']
+    noiseStrength = argdict['noiseStrength']
     
     # Retained for debugging
     isStochastic = True
@@ -324,7 +326,7 @@ def simulateAndSample(argdict):
                                      genelist, proteinlist,
                                      varmapper,revvarmapper)
         
-        P = simulator.simulateModel(Model, y0_exp, pars, isStochastic, tspan, seed)
+        P = simulator.simulateModel(Model, y0_exp, pars, isStochastic, tspan, seed, noiseStrength)
         P = P.T
         retry = False
         ## Extract Time points
diff --git a/BoolODE/simulator.py b/BoolODE/simulator.py
index e993881..d7cd2aa 100644
--- a/BoolODE/simulator.py
+++ b/BoolODE/simulator.py
@@ -1,10 +1,14 @@
 import numpy as np
 
-def noise(x,t):
-    # Controls noise proportional to
-    # square root of activity
-    c = 10.#4.
-    return (c*np.sqrt(abs(x)))
+def noise(x,t, noiseStrength):
+    """
+    Controls noise proportional to square root of activity.
+    default setting of noise strength is 10, defined in JobSettings.
+    NOTE: We have previously tried the value of 4, it might be worth doing 
+    a scan of parameter values to see how sensitive BoolODE output is to 
+    this value.
+    """
+    return (noiseStrength*np.sqrt(abs(x)))
 
 def deltaW(N, m, h,seed=0):
     """Generate sequence of Wiener increments for m independent Wiener
@@ -17,7 +21,7 @@ def deltaW(N, m, h,seed=0):
     np.random.seed(seed)
     return np.random.normal(0.0, h, (N, m))
 
-def eulersde(f,G,y0,tspan,pars,seed=0.,dW=None):
+def eulersde(f,G,y0,tspan,pars,noiseStrength, seed=0., dW=None):
     """
     Adapted from sdeint implementation https://github.com/mattja/sdeint/
 
@@ -53,7 +57,7 @@ def eulersde(f,G,y0,tspan,pars,seed=0.,dW=None):
         tn = currtime
         yn = y[n]
         dWn = dW[n,:]
-        y[n+1] = yn + f(yn, tn,pars)*h + np.multiply(G(yn, tn),dWn)
+        y[n+1] = yn + f(yn, tn,pars)*h + np.multiply(G(yn, tn, noiseStrength),dWn)
         # Ensure positive terms
         for i in range(len(y[n+1])):
             if y[n+1][i] < 0:
@@ -62,7 +66,7 @@ def eulersde(f,G,y0,tspan,pars,seed=0.,dW=None):
         n += 1 
     return y
 
-def simulateModel(Model, y0, parameters,isStochastic, tspan,seed):
+def simulateModel(Model, y0, parameters, isStochastic, tspan, seed, noiseStength):
     """Call numerical integration functions, either odeint() from Scipy,
     or simulator.eulersde() defined in simulator.py. By default, stochastic simulations are
     carried out using simulator.eulersde.
@@ -87,7 +91,7 @@ def simulateModel(Model, y0, parameters,isStochastic, tspan,seed):
     if not isStochastic:
         P = odeint(Model,y0,tspan,args=(parameters,))
     else:
-        P = eulersde(Model,noise,y0,tspan,parameters,seed=seed)
+        P = eulersde(Model,noise,y0,tspan,parameters,noiseStength, seed=seed)
     return(P)
 
 def getInitialCondition(ss, ModelSpec, rnaIndex,
diff --git a/config-files/example-config.yaml b/config-files/example-config.yaml
index 5fc7da7..2c9af64 100644
--- a/config-files/example-config.yaml
+++ b/config-files/example-config.yaml
@@ -107,6 +107,11 @@ jobs:
     ## by a smaller threshold of activation. Specify the strength
     ## as a number between 1 and 20. (Max protein level=20)
     interaction_strengths: "dyn-linear_strengths.txt"
+
+    ################ EXPERIMENTAL FEATURES ################
+    ## Noise strength to use in stochastic simulations
+    # noise_strength: 10
+
     
 post_processing:
   ## Once the simulations have been performed, individual trajectories are