# examples Package¶

## examples Package¶

Examples for the cgptoolbox.

## basic Module¶

Basic example of cGP study of the FitzHugh-Nagumo model of nerve membrane.

from cgp.examples.basic import fitzhugh
import scipy.integrate

t = np.arange(0, 50, 0.1)
y = scipy.integrate.odeint(fitzhugh, [-0.6, -0.6], t)
plt.plot(t, y)
plt.legend(["V (transmembrane potential)", "W (recovery variable)"])


The code in this file is ready to run, but is primarily written to be read. It exemplifies the building blocks of a cGP model study:

Component Description Example
Generation of genotypes Genotypic variation whose phenotypic effects are to be estimated. Full enumeration feasible if the genotype space is small. Other options include pedigrees or statistical experimental designs.
Genotype-to-parameter map Lumping together lower-level physiology such as gene regulation Maximum conductances of ion channels
Parameter-to-phenotype map Physiological model Heart cell electrophysiology
Virtual experiment Manipulating the physiological model Pacing with electrical stimuli
Aggregate phenotypes Summarizing the raw model output Action potential duration
Analysis, summary, visualization

This very simple example defines most components in one file for illustration. Real applications will typically draw upon existing code libraries, model repositories and databases. Linking these together may require a lot of code, a complication which we ignore here to make the main concepts stand out more clearly.

Our example system is the FitzHugh-Nagumo model of an excitable heart muscle cell. It has four parameters, one of which represents a stimulus current. Virtual pacing of the cell can be implemented by setting the stimulus current to a nonzero value at regular intervals. The resulting “action potential” (time-course of transmembrane voltage) is often characterized by its duration, for instance APD90 for the action potential duration to 90 % repolarization.

Here, we assume that each parameter is governed by a single gene of which there are two alleles. This trivial genotype-to-parameter map is a caricature in the absence of actual data, but is a minimal assumption to make the model amenable to cGP analysis.

Here we can compute phenotypes for all 27 genotypes (three loci with three possibilities each gives aa bb cc, aa bb Cc, aa bb CC, aa Bb cc, ..., AA BB CC). In higher-dimensional cases, the genotype space must be sampled according to some experimental design, possibly constrained by pedigree information.

References:

In the original version (FitzHugh 1961), the definition of the transmembrane potential is such that it decreases during depolarization, so that the action potential starts with a downstroke, contrary to the convention used in FitzHugh 1969 and in most other work. The equations are also somewhat rearranged. However, figure 1 of FitzHugh 1961 gives a very good overview of the phase plane of the model.

fitzhugh(y, t=None, par=rec.array([(0.7, 0.8, 0.08, 0.0)],
dtype=[('a', '<f8'), ('b', '<f8'), ('theta', '<f8'), ('I', '<f8')]))

FitzHugh (1969) version of FitzHugh-Nagumo model of nerve membrane.

Default parameter values are from FitzHugh (1969), figure 3-2.

Parameters: y (array_like) – State vector [V, W],where V is the transmembrane potential and W is a “recovery variable”. t (scalar) – Time. Ignored, but required by scipy.integrate.odeint. a, b, theta (float) – Positive constants. I (float) – Transmembrane stimulus current.

References: see module docstring.

cgp.virtexp.examples.Fitz.

cgp.examples.basic.gt2par(gt)[source]

Genotype-to-parameter mapping.

cgp.examples.basic.par2ph(par)[source]

Parameter-to-phenotype mapping.

cgp.examples.basic.ph2agg(ph, tol=0.001)[source]

Phenotype aggregation.

cgp.examples.basic.summarize(gt, agg)[source]

Scatterplot matrix of aggregated phenotypes vs genotypes.

Adjust axis limits to fully show markers.

cgp.examples.basic.cgpstudy()[source]

Basic example of connecting the building blocks of a cGP study.

This top-level orchestration can be done in many different ways, depending on personal preference and on the need for features such as storage, caching, memory management or parallelization. Some examples are given in cgp.examples.hpc.

## hpc Module¶

Demo of infrastructure for large-scale cGP studies.

cgp.examples.hpc.hdfcaching()[source]

Auto-cache/save results to HDF.

cgp.examples.hpc.clusterjob()[source]

Splitting tasks as arrayjobs on a PBS cluster.

## sensitivity Module¶

Sensitivity analysis of action potential duration, bridging R and Python.

class cgp.examples.sensitivity.Model(workspace='bondarenko_szigeti_bett_kim_rasmusson_2004', *args, **kwargs)[source]
cgp.examples.sensitivity.mat2par(mat)[source]

Make parameter recarray from R matrix.

## sensitivity_beeler Module¶

Sensitivity analysis of action potential duration, bridging R and Python.

class cgp.examples.sensitivity_beeler.Model(url=None, workspace=None, exposure=None, changeset=None, variant=None, localfile=None, t=[0, 1], y=None, p=None, rename={}, use_cython=True, purge=False, **kwargs)[source]

Mix and match virtual experiments.

cgp.examples.sensitivity_beeler.mat2par(mat)[source]

Make parameter recarray from R matrix.

cgp.examples.sensitivity_beeler.scalar_pheno(field)[source]

Make a function to return a named field of the phenotype array.

## sigmoid Module¶

Basic example of cGP study using sigmoid model of gene regulatory networks.

This example is taken from Gjuvsland et al. (2011), the gene regulatory model is eq.(15) in that paper, and the example code reproduces figure 4a.

from cgp.examples.sigmoid import cgpstudy, plot_gpmap
gt, par, ph = cgpstudy()
plot_gpmap(gt, ph)


The gene regulatory model is a Sigmoidmodel. It is a diploid model, in the sense that there is two differential equations per gene, one for each allele. This calls for a genotype-to-parameter map where the two alleles are mapped to two independent sets of parameter values. The function geno2par_diploid() offers such genotype-to-parameter mapping with allele-specific parameters .

Sigmoidmodel encodes a system with three genes and the parameter adjmatrix is used to specify connetivity and mode of regulation. The model in Gjuvsland et al. (2011) only has two genes, so we focus only on gene 1 and 2, and make sure that gene 3 is not regulating any of them.

For a given parameter set the resulting phenotype is the equilibrium expression level of gene 2. In order to find the equilibrium we combine Sigmoidmodel with AttractorMixin.

Reference:

Model(t=(0, 1), y=array([ 0., 0., 0., 0., 0., 0.]), p=rec.array([ ([182.04, 159.72], [20.0, 20.0], [5.5, 5.5], [20.0, 20.0], [5.5, 5.5], [10.0, 10.0], [190.58, 135.08], [32.06, 20.86], [8.84, 7.91], [20.0, 20.0], [5.5, 5.5], [10.0, 10.0], [150.0, 150.0], [20.0, 20.0], [5.5, 5.5], [20.0, 20.0], [5.5, 5.5], [10.0, 10.0])],
dtype=[('alpha1', '<f8', (2,)), ('theta11', '<f8', (2,)), ('p11', '<f8', (2,)), ('theta12', '<f8', (2,)), ('p12', '<f8', (2,)), ('gamma1', '<f8', (2,)), ('alpha2', '<f8', (2,)), ('theta21', '<f8', (2,)), ('p21', '<f8', (2,)), ('theta22', '<f8', (2,)), ('p22', '<f8', (2,)), ('gamma2', '<f8', (2,)), ('alpha3', '<f8', (2,)), ('theta31', '<f8', (2,)), ('p31', '<f8', (2,)), ('theta32', '<f8', (2,)), ('p32', '<f8', (2,)), ('gamma3', '<f8', (2,))]), adjmatrix=array([[0, 0, 0],
[1, 0, 0],
[0, 1, 0]]), boolfunc=None)

Bases: cgp.sigmoidmodels.sigmoid_cvode.Sigmoidmodel, cgp.phenotyping.attractor.AttractorMixin

Class that combines a sigmoid model with the ability to find equilibria.

cgp.examples.sigmoid.par2ph(par)[source]

Phenotype: Equilibrium expression level of gene 2.

This maps a parameter value to a phenotype value.

cgp.examples.sigmoid.cgpstudy()[source]

Connecting the building blocks of a cGP study.

This implements the simulation pipeline in Gjuvsland et al. (2011).

cgp.examples.sigmoid.plot_gpmap(gt, ph)[source]

Lineplot of genotypes and corresponding phenotypes.

## simple Module¶

A very simple example of a causally cohesive genotype-phenotype study.

The code in this file is ready to run, but is primarily written to be read. It exemplifies the building blocks of a cGP model study:

Component Description Example
Generation of genotypes Genotypic variation whose phenotypic effects are to be estimated. Full enumeration feasible if the genotype space is small. Other options include pedigrees or statistical experimental designs.
Genotype-to-parameter map Lumping together lower-level physiology such as gene regulation Maximum conductances of ion channels
Parameter-to-phenotype map Physiological model Heart cell electrophysiology
Virtual experiment Manipulating the physiological model Pacing with electrical stimuli
Aggregate phenotypes Summarizing the raw model output Action potential duration
Analysis, summary, visualization

This very simple example defines everything in one file for illustration. Real applications will typically draw upon existing code libraries, model repositories and databases. Linking these together may require a lot of code, a complication which we ignore here to make the main concepts stand out more clearly.

Our example system is the FitzHugh-Nagumo model of an excitable heart muscle cell. It has four parameters, one of which represents a stimulus current. Virtual pacing of the cell can be implemented by setting the stimulus current to a nonzero value at regular intervals. The resulting “action potential” (time-course of transmembrane voltage) is often characterized by its duration, for instance APD90 for the action potential duration to 90 % repolarization.

Here, we assume that each parameter is governed by a single gene of which there are two alleles. This trivial genotype-to-parameter map is a caricature in the absence of actual data, but is a minimal assumption to make the model amenable to cGP analysis.

Here we can compute phenotypes for all 27 genotypes (three loci with three possibilities each gives aa bb cc, aa bb Cc, aa bb CC, aa Bb cc, ..., AA BB CC). In higher-dimensional cases, the genotype space must be sampled according to some experimental design, possibly constrained by pedigree information.

Large-scale computations will benefit from caching and parallelization. Examples are given in simple_joblib, simple_ipython_parallel, simple_hdfcache.

cgp.examples.simple.monogenicpar(genotype, hetpar, relvar=0.5, absvar=None)[source]

Genotype-to-parameter map that assumes one biallelic locus per parameter.

Parameters: genotype – sequence where each item is 0, 1, or 2, denoting the “low” homozygote, the “baseline” heterozygote, and the “high” homozygote, respectively. hetpar (recarray) – parameter values for “fully heterozygous” individual. relvar (float) – proportion change. absvar (array_like) – absolute change, overrides relvar if present Record array of parameter values for the given genotype.

Gene/parameter names are taken from the fieldnames of hetpar.

Example: Define the baseline parameter values for the full heterozygote, as well as an example genotype whose three loci are high homozygous, heterozygous, and low homozygous, respectively. Mapping this genotype onto parameter space, we see that the parameter values are 150%, 100% and 50% of their respective baselines, as assumed.

>>> hetpar = np.rec.fromrecords([(0.75, 0.5, 0.25)],
...     names=["a", "b", "theta"])
>>> genotype = [2, 1, 0]
>>> monogenicpar(genotype, hetpar)
rec.array([(1.125, 0.5, 0.125)],
dtype=[('a', '<f8'), ('b', '<f8'), ('theta', '<f8')])


Specifying relative or absolute parameter variation.

>>> monogenicpar(genotype, hetpar, relvar = 1.0)
rec.array([(1.5, 0.5, 0.0)],...
>>> monogenicpar(genotype, hetpar, absvar = [1.5, 0.5, 0.0])
rec.array([(2.25, 0.5, 0.25)], ...

cgp.examples.simple.fitzhugh(y, _t=None, a=0.7, b=0.8, theta=0.08, I=0.0)[source]

FitzHugh (1969) version of FitzHugh-Nagumo model of nerve membrane

Default parameter values are from FitzHugh (1969), figure 3-2.

Parameters: y (array_like) – State vector [V, W],where V is the transmembrane potential and W is a “recovery variable”. _t (scalar) – Time. Ignored, but required by scipy.integrate.odeint. a, b, theta (float) – Positive constants. I (float) – Transmembrane stimulus current.
from cgp.examples.simple import *
import scipy.integrate

t = np.arange(0, 50, 0.1)
y = scipy.integrate.odeint(fitzhugh, [-0.6, -0.6], t)
plt.plot(t, y)
plt.legend(["V (transmembrane potential)", "W (recovery variable)"])


References:

In the original version (FitzHugh 1961), the definition of the transmembrane potential is such that it decreases during depolarization, so that the action potential starts with a downstroke, contrary to the convention used in FitzHugh 1969 and in most other work. The equations are also somewhat rearranged. However, figure 1 of FitzHugh 1961 gives a very good overview of the phase plane of the model.

cgp.examples.simple.eq(func, Y0, disp=False, *args, **kwargs)[source]

Find equilibrium by minimizing the norm of the rate vector.

Parameters: func (function) – Function to be minimized. Y0 (array_like) – Starting point for estimation. disp (bool) – Print diagnostics from the minimization? **kwargs (*args,) – Additional arguments passed to func().
>>> eq(fitzhugh, [0, 0])
array([-1.19942411, -0.62426308])

cgp.examples.simple.ap(func, Y0, t0=0.0, stim_period=100.0, stim_amplitude=0.7, stim_duration=1.0, dt=1.0, curr_name='I', *args, **kwargs)[source]

Stimulus-induced action potential.

Parameters: func (function) – Right-hand side of ordinary differential equation which can be passed to scipy.integrate.odeint(). Y0 (ndarray) – Initial state vector. t0 (float) – Initial time; helpful in keeping cumulative time for successive action potentials. stim_period (float) – Time between start of successive stimuli. stim_amplitude (float) – Amplitude of stimulus current. stim_duration – Duration of stimulus current. dt (float) – Time resolution of integration. curr_name (str) – Name of the stimulus current argument. **kwargs (*args,) – Passed to func(). (t, y), where y[i] is state at time t[i].
from cgp.examples.simple import *
t, Y = ap(fitzhugh, [-1.2, -0.6])
plt.subplot(121)
plt.plot(t, Y)
plt.legend(["V", "W"])
plt.subplot(122)
plt.plot(*Y.T)
plt.xlabel("V")
plt.ylabel("W")

cgp.examples.simple.aps(func, Y0, n=5, apfunc=<function ap at 0x0D219C70>, *args, **kwargs)[source]

Consecutive stimulus-induced action potentials.

Returns a list with one (time, states) tuple per action potential.

You can iterate over the list of tuples like so:

>>> for t, Y in aps(fitzhugh, [0, 0]): # iterating over list of tuples
...     pass # e.g. pylab.plot(t, Y)


The list of tuples can be unpacked directly if n is known.

>>> (t0, Y0), (t1, Y1) = aps(fitzhugh, [0, 0], n=2)


Time is reckoned consecutively, not restarting for each action potential.

>>> t0[-1] == t1[0]
True


The default parameter values for fitzhugh() give a regular action potential. Shortening the period can give alternans, whereas an insufficient stimulus will not elicit an action potential.

cgp.examples.simple.recovery_time(v, t=None, p=0.9)[source]

Time to p*100% recovery from first local extremum to initial value.

Intended for action potential duration and the like.

Parameters: t (array_like) – Time, defaults to [0, 1, 2, ...].

Examples:

>>> recovery_time([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 5, 0])
11.8...
>>> recovery_time([10, 0, 20], [100, 101, 102])
101.45
>>> recovery_time([10, 0, 20], p=0.5)
1.25

cgp.examples.simple.ph(gt)[source]

Ad hoc function to compute the phenotype of a single genotype

cgp.examples.simple.visualize(result)[source]

Trivial example of visualization.

## simple_hdfcache Module¶

Cache computations to HDF file for simple genotype-phenotype example.

Here, the line “ph = hdfcache.cache(ph)” is equivalent to the decorator syntax:

@hdfcache.cache
def ph(...):
...


However, cgp.examples.simple.ph is left undecorated so we can illustrate different tools for caching and parallelization.

## simple_ipython_parallel Module¶

Using IPython.parallel to parallelize the computations in simple.py.

ipcluster start must be run for this example to work.

Parallelization can be conveniently applied using decorator syntax:

@lv.parallel()
def f(x):
...


However, cgp.examples.simple.ph is left undecorated so we can illustrate different tools for caching and parallelization.

## simple_joblib Module¶

Using joblib to cache and parallelize the computations in simple.py.

Caching can be conveniently applied using decorator syntax:

@mem.cache
def f(x):
...


However, parallelization won’t work with this syntax as of joblib 0.6.3.