fcmaes - a Python 3 gradient-free optimization library

Parameter Sweeping for Biochemical Reactions

This tutorial

Is related to VilarOscillator.py and Large-scale simulation-based experiments with stochastic models using machine learning assisted approaches.
Explains semi-automated parameter sweeping for stochastic biochemical reaction networks.
Shows how to use multi-objective optimization to:
- Parallelize complex expensive stochastic simulations.
- Guide the parameter sweep toward specific network properties.

The code for this tutorial is here: vilar.py

Motivation

Stochastic simulations of biochemical reaction networks can capture intrinsic noise. This matters when the number of molecular species is low. These simulations are expensive, and the noise makes parameter evaluation harder. For this reason, Fredrik Wrede proposes a machine-learning-based interactive method for model exploration through parameter sweeps on a multi-node CPU network. See the corresponding video.

Our idea is similar, but it uses a different kind of human interaction. Here the user designs a multi-objective fitness function, and optimization then guides the parameter sweep. The goal is not to find a single optimum. Instead, we introduce some "bias" into the sweep so it moves toward "interesting" model properties. This is especially useful in high-dimensional parameter spaces, where exhaustive scanning is not practical. The two approaches could also be combined. The time series generated during optimization could feed the machine-learning-based interactive method.

The Vilar Oscillator

In Mechanisms of noise-resistance in genetic oscillators Jose M.G.Vilar presented a biochemical model of a "circadian clock" that lets organisms maintain an internal sense of daily time. This model can be simulated using GillesPy2, see VilarOscillator.py. The Vilar model has 15 parameters, and we want to understand:

Is the oscillating behavior of the model dependent on specific parameter settings?
Can we find parameters which can affect the oscillating property of the model negatively?
Or does the model have "self-regulating" properties preserving the steady oscillation?

Multi-Objective Function

The core idea is to use multi-objective optimization, which requires a fitness function. The question is how to analyze the result of one simulation run, which produces time series for each species, so that we can characterize the "steady oscillation" property.

    # multi processing result list
    results = mp.Manager().list()

    class fcmaes_problem():

        def __init__(self):
            self.bounds = get_bounds(VilarOscillator(), 100)
            self.dim = len(self.bounds.ub)

        def fitness(self, x):
            model = VilarOscillator()
            set_params(model, x)
            res = model.run(algorithm = "SSA")
            # store params, result tuple
            results.append((x, res))
            R = res['R'] # time series for R
            r_mean = np.mean(R)
            r_over = np.array([r for r in R if r > r_mean])
            ilocs_max = argrelextrema(r_over, np.greater_equal, order=3)[0]
            freq = len(ilocs_max) / len(R)
            peak_dists = np.array([ilocs_max[i] - ilocs_max[i-1] for i in range(1, len(ilocs_max))])
            sdev_peak_dist = np.std(peak_dists)
            peaks = (r_over - r_mean)[ilocs_max]
            sdev_amp = np.std(peaks)
            # maximize sdev_peak_dist and sdev_amp
            return [-sdev_peak_dist, -sdev_amp, freq]----

There are more sophisticated options, for example FFT-based analysis. Here we keep it simple and use scipys argrelextrema to identify maxima in the R species. From these peaks we compute the standard deviation of the amplitude and the standard deviation of the peak time distances. Small values for these standard deviations indicate a steady oscillation, so we use them as objectives to maximize. We also add "frequency" as a third objective. This increases the set of non-dominated solutions, the so-called pareto front.

The optimizer can handle different scalings of both parameters and objectives. We could also define objectives as constraints to exclude specific model properties from the sweep. Constraints are given priority while they are violated and are ignored once they are satisfied. Calling a multi-objective optimizer and plotting the result is straightforward:

    popsize = 64 # population size of the evolutionary algorithm
    max_evaluations = popsize*16 # maximum number of evaluation
    # popsize = 256 # population size of the evolutionary algorithm
    # max_evaluations = popsize*96 # maximum number of evaluation
    nobj = 3 # number of objectives
    ncon = 0 # number of constraints
    # stores all values; if capacity is reached, content is replaced by the pareto front
    store = mode.store(problem.dim, nobj, capacity=max_evaluations)
    # perform the parameter sweep by multi objective optimization
    xs, ys = mode.minimize(mode.wrapper(problem.fitness, 3, interval=64, store=store),
                                    nobj, ncon,
                                    problem.bounds, popsize = popsize,
                                    max_evaluations = max_evaluations,
                                    nsga_update=True, workers=min(popsize, mp.cpu_count()))
    # save last population of the evolutionary algorithm
    np.savez_compressed("sweep", xs=xs, ys=ys)
    # save all evaluations
    xs, ys = store.get_xs(), store.get_ys()
    np.savez_compressed("sweep_all", xs=xs, ys=ys)
    # show results
    for x, res in results[:10]:
        print(list(x), list(res['R']))
    moretry.plot("sweep_all", ncon, xs, ys) # plot 2d
    plot3d(xs, ys, "sweep_3d") # plot 3d

The workers parameter enables parallel simulations through Python multiprocessing. On a 16 core AMD 5950 CPU, this gives about 8 simulations per second. This example stays on one machine, but the same approach could also use dask to distribute optimization runs across multiple nodes, as in this paper, and then merge the results.

The store object keeps all parameter settings together with their objective values. We set capacity=max_evaluations so every evaluation fits. store resolves parallel access conflicts with a multiprocessing.Lock. The simulation time series and their corresponding parameter settings are collected in results = mp.Manager().list().

The computed pareto front shows that the maximal standard deviation of the amplitude is around 260, while the maximal standard deviation of the peak time distance is only around 2.95:

[-2.9529177407085023, -95.4039656757866, 0.04488778054862843] [3963.55534, 5.0, 0.07189, 291.00951, 2333.53036, 176.64492, 100.0, 769.04945, 45.31349, 11.30336, 20.92418, 39.33244, 100.0, 3512.29386, 3351.37915]
[-2.951745715363383, -105.61892834474017, 0.04488778054862843] [3192.0246, 21283.45224, 0.40744, 2282.94153, 2449.13288, 14.40633, 85.15188, 1000.0, 8.95256, 7.73864, 71.91611, 0.02, 10.65094, 908.59928, 7486.3519]
[-2.9056680082890516, -118.76987812832546, 0.04488778054862843] [4750.8944, 21187.26512, 0.30343, 3053.65532, 958.32784, 500.0, 2.02548, 0.1, 17.85172, 5.61477, 70.8133, 101.64201, 38.998, 0.5, 2709.01493]
[-2.829650228424869, -126.14233583095285, 0.04488778054862843] [1844.98312, 8018.1519, 0.39544, 1950.09358, 199.70482, 0.05, 22.21654, 776.71529, 19.2612, 14.99067, 43.73244, 146.38544, 61.3209, 2707.33989, 9282.14773]
[-2.8272034919289215, -133.07568857191342, 0.04488778054862843] [2727.76445, 50000.0, 0.50671, 1304.37028, 2589.98677, 464.17394, 85.61988, 630.31469, 35.83535, 4.50604, 89.76378, 49.71666, 39.89323, 975.28531, 10000.0]
[-2.6299556396765835, -135.47628100402463, 0.04738154613466334] [5000.0, 30524.08898, 0.73944, 3784.11623, 252.9115, 500.0, 37.03655, 832.33873, 12.4825, 7.32541, 39.97521, 78.18922, 82.90321, 4001.87626, 7163.73231]
[-2.55255862921022, -101.80509488110111, 0.0399002493765586] [1679.00163, 49172.18683, 0.0001, 698.72041, 189.6883, 231.44121, 91.28288, 827.96068, 19.55392, 9.63559, 59.86862, 57.8062, 75.01048, 5000.0, 2953.52827]
[-2.4462982238476156, -157.9285458855882, 0.04239401496259352] [1616.84502, 25773.2103, 0.19491, 3311.77833, 1702.56557, 72.51257, 48.11612, 297.08423, 18.25544, 11.40069, 31.05834, 0.02, 18.77209, 2323.45674, 2993.41434]
[-2.434930982366623, -116.0010607440725, 0.0399002493765586] [2909.06717, 8230.83586, 0.34597, 4137.00498, 0.5, 0.05, 53.08472, 0.1, 0.005, 16.94276, 28.20045, 45.02235, 22.52807, 127.21803, 1316.2461]
[-2.3795424396766327, -116.71493049306073, 0.0399002493765586] [0.5, 38009.30647, 0.82832, 4025.69335, 2842.63817, 426.54968, 29.56595, 284.86072, 38.63421, 0.002, 31.32197, 200.0, 95.21752, 0.5, 10000.0]
[-2.3570226039551585, -168.11016430825947, 0.0399002493765586] [2226.31818, 50000.0, 1.0, 0.5, 2520.63258, 63.26119, 91.53521, 131.19077, 32.69053, 6.16228, 100.0, 89.46733, 61.56866, 741.69312, 6533.65796]
[-2.314390067400706, -207.69080176021873, 0.04488778054862843] [1738.84673, 28325.49892, 0.0001, 4153.941, 0.5, 297.23208, 47.08712, 285.88175, 18.25992, 0.002, 75.07047, 127.18269, 40.57461, 0.5, 9832.98613]
[-2.1758618981911515, -180.4157060502754, 0.04239401496259352] [5000.0, 50000.0, 0.09536, 1102.04172, 507.11942, 131.34008, 23.0144, 825.01715, 22.50374, 10.67058, 0.01, 0.02, 92.29957, 1366.8583, 965.63167]
[-1.3942300925673667, -114.8690849039317, 0.03740648379052369] [3966.74083, 8256.54672, 1.0, 827.22444, 3410.66379, 482.81114, 23.71309, 643.26172, 28.77231, 20.0, 42.67218, 0.02, 88.48254, 4560.1138, 7150.24457]
[-1.3420210162897097, -146.56429760802368, 0.03740648379052369] [842.00082, 50000.0, 0.19752, 650.88578, 1492.86002, 64.03229, 100.0, 630.28532, 15.78265, 9.23157, 0.01, 25.39473, 58.1836, 1795.3776, 3223.63713]
[-1.3228756555322954, -216.46190839375674, 0.04239401496259352] [3316.9303, 49172.18683, 0.59904, 698.72041, 189.6883, 0.05, 100.0, 827.96068, 19.55392, 9.63559, 59.86862, 153.03553, 97.46121, 5000.0, 2953.52827]
[-1.2892719737209144, -186.11068615208526, 0.0399002493765586] [2380.31009, 36173.28215, 0.0001, 4032.86154, 2955.06148, 300.79915, 76.97691, 976.55334, 50.0, 0.002, 65.68229, 100.41598, 53.51887, 4936.06836, 7099.12308]
[-1.263700429496724, -147.90483426852552, 0.03740648379052369] [5000.0, 41614.45492, 0.0001, 2042.40438, 3849.19248, 155.90764, 58.34843, 364.50473, 13.89905, 2.23895, 100.0, 152.30698, 55.85457, 4662.94334, 1.0]
[-1.254325848148452, -190.04792775442198, 0.0399002493765586] [3981.57454, 40087.32536, 0.307, 319.81802, 2846.01873, 188.623, 40.72906, 522.07154, 0.005, 11.26821, 35.0033, 97.53761, 19.81748, 2712.01541, 6838.85863]
[-1.1925695879998877, -198.70734605192632, 0.0399002493765586] [5000.0, 25412.65676, 0.0001, 3679.39956, 2835.75002, 208.42794, 97.36303, 894.59821, 23.4597, 4.90817, 19.93437, 95.33524, 15.36248, 2324.47721, 10000.0]
[-1.178030178747903, -151.3895049936495, 0.03740648379052369] [1363.62001, 31011.23828, 0.48028, 4863.33951, 963.52964, 213.70395, 30.89808, 358.17736, 50.0, 13.69323, 57.41104, 109.74502, 88.74388, 1514.7859, 5920.81353]
[-1.0873004286866728, -211.57766747698116, 0.0399002493765586] [0.5, 44679.14821, 0.65303, 1910.94801, 2314.8589, 173.2861, 55.17358, 367.67874, 0.005, 1.1258, 51.20502, 157.59849, 100.0, 4478.88327, 1.0]
[-1.0749676997731399, -222.04714272818285, 0.0399002493765586] [3247.96439, 1352.80086, 0.13737, 2351.74198, 2294.67424, 279.65285, 98.58432, 503.23886, 26.17302, 15.31982, 47.78622, 87.54788, 100.0, 2344.1427, 2758.37033]
[-1.0326308782000686, -194.84854517177067, 0.03740648379052369] [459.55155, 50000.0, 0.76135, 4863.33951, 963.52964, 213.70395, 24.98266, 242.57983, 50.0, 2.9316, 57.41104, 200.0, 96.65515, 1514.7859, 10000.0]
[-1.019803902718557, -234.11829274962687, 0.0399002493765586] [5000.0, 44117.05737, 0.90691, 4137.00498, 0.5, 138.02771, 98.00668, 346.56105, 13.56542, 17.50474, 34.12235, 0.02, 19.33328, 492.70916, 10000.0]
[-0.9354143466934853, -236.54461820422816, 0.04239401496259352] [3932.05576, 35815.82508, 1.0, 1955.83674, 2910.72433, 169.00539, 11.68632, 0.1, 50.0, 12.20671, 98.82335, 0.02, 67.57933, 3636.91003, 9.92914]
[-0.9333240326917549, -254.47453152964408, 0.04239401496259352] [2814.62758, 50000.0, 0.58636, 3204.78046, 1270.30639, 135.21917, 29.06001, 917.14848, 39.93926, 0.002, 0.01, 104.66887, 73.35097, 382.12618, 3387.31236]
[-0.9285592184789413, -244.861113376032, 0.0399002493765586] [3670.84244, 822.05728, 0.201, 3462.9566, 1997.32353, 256.3988, 86.54088, 875.1952, 0.005, 4.11743, 34.55674, 84.21904, 37.88851, 2254.81908, 1.0]
[-0.8637671850678283, -256.2949884552299, 0.04239401496259352] [4185.77373, 22788.5206, 0.53141, 0.5, 0.5, 173.65096, 46.7135, 252.04982, 24.21218, 13.53778, 82.29212, 104.30263, 82.92906, 887.3639, 7338.28102]
[-0.618241233033047, -258.4637305309973, 0.0399002493765586] [2860.71177, 29392.31481, 0.55805, 3573.22623, 0.5, 357.77082, 85.8377, 515.86832, 45.87822, 6.8493, 48.78439, 90.04616, 70.66459, 245.69653, 10000.0]

In particular, the standard deviation of the peak time distances stays stable, although we do see some amplitude deviations.

Now let us try the opposite direction. We want to identify parameters that stabilize the oscillation and minimize both standard deviations:

        def fitness(self, x):
            ...
            # minimize sdev_peak_dist and sdev_amp
            return [sdev_peak_dist, sdev_amp, freq]

This time the parameter sweep is guided in the opposite direction, and it finds many low-standard-deviation solutions:

[0.33993463423951903, 85.05106462443608, 0.0399002493765586] [4667.13481, 46359.99664, 0.50975, 4963.54422, 11.66221, 0.05, 51.12464, 402.27105, 10.53585, 1.43816, 27.0709, 133.96803, 100.0, 2996.65492, 7754.62117]
[0.4, 74.98656129600823, 0.0399002493765586] [4085.55047, 41534.06242, 0.25192, 266.76638, 4828.90463, 421.50792, 59.04809, 0.1, 26.36669, 16.87143, 60.35491, 86.88787, 78.07726, 2996.71963, 4881.88487]
[0.41032590332414487, 114.08737392406263, 0.03740648379052369] [5000.0, 50000.0, 0.23102, 0.5, 0.5, 249.3216, 0.01, 889.09787, 50.0, 20.0, 50.98781, 149.82714, 0.01, 3304.80047, 38.29661]
[0.410325903324145, 90.36125029876222, 0.03740648379052369] [3922.44784, 50000.0, 0.47131, 5000.0, 225.90615, 268.70649, 59.34419, 955.78435, 50.0, 20.0, 56.17379, 52.21683, 0.01, 2545.12168, 7786.30862]
[0.4422166387140533, 70.73870824202262, 0.0399002493765586] [0.5, 8884.6834, 0.14039, 2121.80486, 5000.0, 9.94685, 0.01, 515.15081, 31.46132, 9.10421, 81.49082, 137.09595, 0.01, 397.43259, 1.0]
[0.45175395145262565, 89.4535012667972, 0.03740648379052369] [3969.552, 35595.20059, 0.47131, 5000.0, 225.90615, 500.0, 59.34419, 669.16249, 45.03994, 4.21045, 56.17379, 52.21683, 18.96992, 522.35537, 8620.9769]
[0.4573660169594892, 83.93870779721752, 0.03740648379052369] [3750.125, 5.0, 0.31521, 1321.85626, 1259.17944, 311.64924, 9.56078, 355.97666, 33.78151, 11.03842, 0.01, 145.53896, 71.08417, 3806.81813, 8101.10279]
[0.47140452079103173, 53.79616476803899, 0.0399002493765586] [4877.2711, 50000.0, 0.37939, 509.91588, 5000.0, 95.47948, 0.01, 1000.0, 50.0, 17.43156, 71.17818, 169.94934, 15.65189, 3889.87649, 38.29661]
[0.4791574237499549, 70.6209285945431, 0.03740648379052369] [2148.41722, 3823.11184, 0.76019, 0.5, 2526.24176, 85.91607, 0.01, 609.97549, 44.41766, 20.0, 65.27637, 49.73482, 26.75771, 73.16505, 10000.0]
[0.48989794855663554, 46.94261789674709, 0.0399002493765586] [4486.38357, 41043.38753, 0.30004, 1153.50072, 1841.51075, 0.05, 48.91541, 388.84663, 34.56827, 19.82925, 46.43581, 194.80536, 11.18324, 857.99832, 3642.3788]
[0.4948716593053934, 65.56794609834562, 0.03740648379052369] [3026.12731, 5.0, 1.0, 4955.40119, 3473.91023, 16.57025, 88.00157, 478.59413, 23.4971, 0.002, 15.32885, 91.05324, 79.769, 315.13728, 10000.0]
[0.498887651569859, 41.541395017500314, 0.0399002493765586] [2401.0277, 7612.15779, 0.31682, 79.56769, 0.5, 471.41371, 95.95611, 336.75528, 3.92945, 17.79122, 7.11661, 191.06115, 36.05853, 299.63769, 152.02355]
[0.5, 63.90305157032801, 0.03740648379052369] [3521.04674, 50000.0, 0.33635, 4672.30343, 4993.83557, 321.52988, 53.77845, 969.06066, 50.0, 20.0, 100.0, 103.33393, 95.49405, 2637.95239, 10000.0]
[0.573488351136175, 41.42293295989553, 0.0399002493765586] [2772.0545, 8312.37421, 0.35097, 3282.00071, 1731.19102, 19.8691, 89.37192, 596.8713, 12.83699, 20.0, 100.0, 200.0, 13.03374, 2080.06046, 3750.0288]
[0.6110100926607787, 40.37712842687058, 0.0399002493765586] [2700.9404, 43773.06096, 0.61493, 93.38473, 4686.28757, 67.37642, 61.25017, 720.21311, 1.1664, 20.0, 36.81052, 171.52996, 39.75133, 1951.84428, 5994.53237]
[0.6226998490772391, 56.23616866987532, 0.03740648379052369] [441.53957, 28266.85623, 0.4996, 0.5, 4241.41337, 321.28582, 50.78861, 1000.0, 9.58572, 12.72225, 63.05839, 36.69093, 88.46005, 4363.74585, 1.0]
[0.6388765649999398, 50.01892974997899, 0.03740648379052369] [2786.56457, 1752.71181, 0.7316, 4583.21976, 0.5, 73.88818, 67.44368, 454.1061, 1.39891, 20.0, 100.0, 200.0, 14.71424, 2165.23854, 9916.8418]
[0.7034898429854359, 44.120239749524885, 0.03740648379052369] [775.25653, 38721.18599, 0.91891, 4604.08801, 4744.04525, 353.08363, 100.0, 901.34442, 16.91242, 2.915, 49.48101, 64.53518, 48.74173, 584.5451, 1388.3146]
[0.7180219742846006, 38.740674281044726, 0.0399002493765586] [1089.1402, 25002.5, 0.65476, 2623.74255, 4494.51243, 187.9306, 21.40189, 991.73702, 2.32626, 6.64121, 96.52458, 79.65669, 60.70325, 3015.74996, 10000.0]
[0.9970370305242863, 114.33385146944336, 0.034912718204488775] [2533.58606, 11377.76193, 0.95959, 245.79493, 3822.78805, 438.11378, 53.71136, 667.52735, 12.34801, 11.26597, 69.41253, 0.02, 43.43411, 595.77608, 7057.76595]
[1.071414482860317, 81.66344572296065, 0.034912718204488775] [5000.0, 5.0, 0.91981, 887.67146, 2719.74535, 286.77574, 11.60726, 0.1, 23.14101, 20.0, 92.87981, 120.10005, 45.65844, 2489.06925, 10000.0]
[1.0986812966989, 79.28378376028063, 0.034912718204488775] [1593.45681, 41602.71265, 0.0001, 1503.00018, 5000.0, 0.05, 27.25116, 504.13417, 50.0, 20.0, 58.40439, 0.02, 9.11069, 3414.72936, 4293.5848]

Visualizing all evaluations

If we visualize all stored evaluations from the optimization-driven parameter sweep, there is not much difference between the two runs. First, the run that maximizes the standard deviations:

Then the run that minimizes the standard deviations:

Finally, a 3D view:

Conclusion

Multi-objective optimization can speed up the parameter sweep of a stochastic biochemical reaction network model.
Simulations are executed in parallel utilizing all processor cores.
The objective function guides the parameter sweep toward "interesting" model properties or, in our example, tries to destroy them.