## Selection of the week

August 15, 2014

Eiene Kleine Nachtmusik by Wolfgang Amadeus Mozart performed by the Concertgebouw Kamerorkest in 2013 in Amsterdam.

## Numerical integration

August 14, 2014

Solving differential equations numerically is a very mature field with multiple algorithms carefully explained in multiple textbooks. However, when it comes down to actually solving them in practice by nonspecialists, a lot of people, myself included, will often resort to the Euler method. There are two main reasons for this. The first is that it is trivial to remember and code and the second is that Moore’s law has increased the speed of computing sufficiently that we can get away with it. However, I think it’s time to get out of the eighteenth century and at least move into the nineteenth.

Suppose you have a differential equation $\dot{x}=f(x,t)$, where $x$ can be a scalar or a vector in an arbitrary number of dimensions.  The Euler method simply discretizes the time derivative via

$\frac{dx}{dt} \rightarrow \frac{x(t+h)-x(h)}{h}$

so that $x$ at time step $t+h$ is given by

$x(t+h)=x(t)+h*f(x(t),t)$

Now what could be simpler than that? There are two problems with this algorithm. The first is that it is only accurate to order $h$ so you need to have very small steps to simulate accurately. The second is that the Euler method is unstable if the step size is too large. The way around these limitations is to use very small time steps, which is computationally expensive.

The second most popular algorithm is probably (4th order) Runge-Kutta, which can generally use larger time steps for the same accuracy. However, a higher order algorithm is not necessarily faster because what really slows down a numerical scheme is the number of function evaluations it needs to make and 4th order Runge-Kutta needs to make a lot of function evaluations. Runge-Kutta is also plagued by the same stability problems as Euler and that can also limit the maximum allowed step size.

Now, these problems have been known for decades if not centuries and work-arounds do exist. In fact, the most sophisticated solvers such as CVODE tend to use adaptive multi-step methods, like Adams-Bashforth (AB) or it’s implicit cousin Adams-Moulton and this is what we should all be implementing in our own codes. Below I will give simple pseudo-code to show how a 2nd order Adams-Bashforth scheme is not much more complicated to code than Euler and uses the same number of function evaluations. In short, with not much effort, you can speed your code immensely by simply switching to AB.

The cleverness of AB is that it uses previous steps to refine the estimate for the next step. In this way, you get a higher order scheme without having to re-evaluate the right hand side of your ODE. You simply reuse your computations. It’s the ultimate green scheme. For the same ODE example above 2nd order AB is simply

$x(t+2h)=x(t+h)+h\left(\frac{3}{2}f(x(t+h),t+h)-\frac{1}{2}f(x(t),t)\right)$

You just need to store the previous two evaluations to compute the next step.  To start the algorithm, you can just use an Euler algorithm to get the first step.

Here is the 2nd order AB algorithm written in pseudo Julia code for a one dimensional version of the ODE above.

# make x a two dimensional vector for the two time steps

# Set initial condition

x[1]=0;

# Store evaluation of function f(x,t) for later use

fstore[1]= f(x[1],h)

# Take Euler step

x[2] = x[1] + h*fstore[1]

# store function at time = 2h

fstore[2]  = f(x[2],2*h)

# Store computed values for output

xout[1]=x[1]

xout[2]=x[2]

# Take AB steps

for t = 3:Tfinal

# set update indices for x based on parity of t

index2 = t%2+1

index1 = (t-1)%2+1

# 2nd order AB step

x[index1] = x[index2]+h*(1.5*fstore[2]-0.5*fstore[1]))

# update stored function

fstore[index1] = f(x[index1],t*h)

xout[t]=x[index1]

end

return xout

## Selection of the Week

August 8, 2014

The 4th movement of Franz Schubert’s Trout Quintet (Piano Quintet in A Major) performed by the Schubert Ensemble.

## The neuro software gap

August 4, 2014

Our foray into Julia was mostly motivated by a desire to have a programming language as fast as C but as easy to program as Matlab. In particular, we need really fast looping speed, which is necessary to simulate dynamical systems. Our particular project involves fitting a conductance-based model (i.e. Hodgkin-Huxley-like equations) to some data using MCMC. We’re not trying to fit the precise spike shapes per se but the firing rates of synaptically coupled pools of neurons in response to various inputs. We thus need to run simulations for some time, compare the averaged rates to previous runs, and iterate. For good fits, we may need to do this a million times. Hence, we need speed and we can’t really parallelize because we need information from one iteration before going to the next.

Now, one would think that simulating neural networks is a long solved problem and I thought that too. Hence, I instructed my fellow Wally Xie to look into all the available neural simulators and see what fits the bill. We first tried Neuron but couldn’t get all of it to work on our Linux machine running CentOS (which is what NIH forces us to use).  That was the end of Neuron for us. We next tried Nest but found it difficult to specify the neuron model. We finally settled on Brian, which is adaptable, written in Python, and supported by very nice people, but is turning out to be too slow. Native Python is also too slow. Stan would be great but it does not yet support differential equations. Wally is now implementing the code in Julia, which is turning out to be very fast but has a possible bug that is preventing long runs. The developers are extremely helpful though so we believe this will be worked out soon. However, it is hard to get plotting to work on Julia and as pointed out by many, Julia is not nearly as complete as Python or Matlab. (Yes, Matlab is too slow). We could use C and will go to it next if Julia does not work out. Thus, even though people have been developing neural simulators for decades, there is still an unmet need for something adaptable, fast, easy to use, and runs out of the box. A neuro-Stan would be great.

## Selection of the Week

August 1, 2014

Summer from Antonio Vivaldi’s Four Seasons of course!  Here is a performance by the 2nd place winner at the 2012 Menuhin Junior violin competition.

## JSMB/SMB talk

July 29, 2014

I just gave an invited plenary talk at the joint annual meeting of the Japanese Society of Mathematical Biology and Society of Mathematical Biology in Osaka Japan.  I talked about my work on steroid-regulated gene transcription.  The slides are here.  Previous posts on the topic, including background summaries, can be found here.

## Tim’s Vermeer

July 28, 2014

Jan Vermeer has been one of my favourite painters ever since I saw his famous “The Girl with the Pearl Earring” painting that was on display in Toronto in the 1980’s. I’ve been on a quest to see all of his paintings although its been on hiatus for the past ten years. Here is the list of what I’ve seen so far (I have five left). You only need to stand in front of a Vermeer for a few seconds to be mesmerized. I stood in front of “The Music Lesson” in Buckingham Palace for at least an hour. The guard started joking with me because I was so transfixed. This is why I’ve been intrigued by recent suggestions by artist David Hockney and others that some great old masters like Vermeer and van Eyck may have used optical aids like camera obscura. Well, inventor Tim Jenison has taken this theory to another level by attempting to completely recreate Vermeer’s Music Lesson using a set up of mirrors and lenses that he (re)invented. The endeavor is documented in the film Tim’s Vermeer directed by Teller of Penn and Teller fame. Whether you believe the theory or not (I actually do and it doesn’t detract at all for my love of Vermeer), what this film does do so well is to show what dedication, thought, patience, and careful execution can accomplish. I got tired just watching him paint the threads in a Persian rug using his optical tool.

## Selection of the week

July 23, 2014

What is called classical music mostly refers to the Western symphony orchestra tradition that starts in the seventeenth century with Vivaldi and peaks in the early twentieth with Mahler. While classical music remains popular, my unscientific sampling of concert hall audiences indicates that the demographic skews to retirement age and above. I don’t know if this means that a generation of music lovers is about to depart or that people only have the patience to sit through a long concert when they are older. In an attempt to introduce a new generation to classical music, I thought I would present a selection each week. And what’s a better way to kick it off then with the pseudo-Baroque precursor to heavy metal, the Praeludium and Allegro by Fritz Kreisler. Kreisler performed in the first half of the twentieth century. He was one of the greatest violin virtuosos of all time and also wrote some great violin ditties. Here is a performance by the then 13 year old Canadian/American violinist Leila Josefowicz in 1991.

## Talk at Jackfest

July 16, 2014

I’m currently in Banff, Alberta for a Festschrift for Jack Cowan (webpage here). Jack is one of the founders of theoretical neuroscience and has infused many important ideas into the field. The Wilson-Cowan equations that he and Hugh Wilson developed in the early seventies form a foundation for both modeling neural systems and machine learning. My talk will summarize my work on deriving “generalized Wilson-Cowan equations” that include both neural activity and correlations. The slides can be found here. References and a summary of the work can be found here. All videos of the talks can be found here.

Addendum: 17:44. Some typos in the talk were fixed.

Addendum: 18:25. I just realized I said something silly in my talk.  The Legendre transform is an involution because the transform of the transform is the inverse. I said something completely inane instead.

## Zero Matlab 2015

July 7, 2014

So I have bitten the bullet and am committed to phasing out Matlab completely by 2015.  I have Julia installed and can sort of run code in it although I have no idea how to change directories.  cd() doesn’t seem to work for me.  I have also tried to install matplotlib on my MacBook Pro running OS 10.8 using a SciPy superpack but it does not seem to work yet.  When I try to plot something, nothing happens but I have some bouncing rocket ships in my dock.  Feel free to let me know what I’m doing wrong.

Addendum: I installed iPython using Pip and I can plot out of iPython

## The need for speed

July 1, 2014

Thanks for all the comments about the attributes of Python and Julia. It seems to me that the most prudent choice is to learn Python and Julia. However, what I would really like to know is just how fast these languages really are and here is the test. What I want to do is to fit networks of coupled ODEs (and  PDEs) to data using MCMC (see here). This means I need a language that loops fast. An example in pseudo-Matlab code would be

for n = 1:N

for i = 1:T

y(i+1) = M\y(i)

end

Compare to data and set new parameters

end

where h is a parameter and M is some matrix (say 1000 dimensional), which is sometimes a Toeplitz matrix but not always. Hence, in each time step I need to invert a matrix, which can depend on time so I can’t always precompute, and do a matrix multiplication. Then in each parameter setting step I need to sum an objective function like the mean square error over all the data points. The code to do this in C or Fortran can be pretty complicated because you have to keep track of all the indices and call linear algebra libraries. I thus want something that has the simple syntax of Matlab but is as fast as C. Python seems to be too slow for our needs but maybe we haven’t optimized the code. Julia seems like the perfect fit but let me know if I am just deluded.

## Julia vs Python

June 29, 2014

I was about to start my trek up Python mountain until Bard Ermentrout tipped me to the Julia language and I saw this speed table from here (lower is faster):

Fortran Julia Python R Matlab Octave Mathe-matica JavaScript Go fib gcc 4.8.1 0.2 2.7.3 3.0.2 R2012a 3.6.4 8.0 V8 3.7.12.22 go1 0.26 0.91 30.37 411.36 1992.00 3211.81 64.46 2.18 1.03 5.03 1.60 13.95 59.40 1463.16 7109.85 29.54 2.43 4.79 1.11 1.14 31.98 524.29 101.84 1132.04 35.74 3.51 1.25 0.86 0.85 14.19 106.97 64.58 316.95 6.07 3.49 2.36 0.80 1.00 16.33 15.42 1.29 237.41 1.32 0.84 1.41 0.64 1.66 13.52 10.84 6.61 14.98 4.52 3.28 8.12 0.96 1.01 3.41 3.98 1.10 3.41 1.16 14.60 8.51

Julia is a dynamic high level language like MATLAB and Python that is open source and developed at MIT. The syntax looks fairly simple and it is about as fast as C (Fortran looks like it still is the Ferrari of scientific computing). Matlab is fast for vector and matrix operations but deadly slow for loops. I had no idea that Mathematica was so fast. Although Julia is still relatively new and not nearly as expansive as Python, should I drop Python for Julia?

June 26, 2014

## New Papers

June 16, 2014
Two new papers are now in print:
The first is on applying compressed sensing to genomics is now published in Gigascience. The summary of the paper is here and the link is here.
The second is on steroid-regulated gene induction and can be obtained here.
Biochemistry. 2014 Mar 25;53(11):1753-67. doi: 10.1021/bi5000178. Epub 2014 Mar 11.

### Abstract

A gene induction competition assay has recently uncovered new inhibitory activities of two transcriptional cofactors, NELF-A and NELF-B, in glucocorticoid-regulated transactivation. NELF-A and -B are also components of the NELF complex, which participates in RNA polymerase II pausing shortly after the initiation of gene transcription. We therefore asked if cofactors (Cdk9 and ELL) best known to affect paused polymerase could reverse the effects of NELF-A and -B. Unexpectedly, Cdk9 and ELL augmented, rather than prevented, the effects of NELF-A and -B. Furthermore, Cdk9 actions are not blocked either by Ckd9 inhibitors (DRB or flavopiridol) or by two Cdk9 mutants defective in kinase activity. The mode and site of action of NELF-A and -B mutants with an altered NELF domain are similarly affected by wild-type and kinase-dead Cdk9. We conclude that Cdk9 is a new modulator of GR action, that Ckd9 and ELL have novel activities in GR-regulated gene expression, that NELF-A and -B can act separately from the NELF complex, and that Cdk9 possesses activities that are independent of Cdk9 kinase activity. Finally, the competition assay has succeeded in ordering the site of action of several cofactors of GR transactivation. Extension of this methodology should be helpful in determining the site and mode of action of numerous additional cofactors and in reducing unwanted side effects.

PMID: 24559102 [PubMed - indexed for MEDLINE]
PMCID: PMC3985961 [Available on 2015/2/21]

## Slides for talk on gene expression

June 10, 2014

Here are the slides for my talk today at the NIH Systems Biology Forum on gene expression.  Background for the talk can be found here.

## Marc Andreesen on EconTalk

June 3, 2014

If you have any interest in technology and the internet then you should definitely listen to this EconTalk podcast with Marc Andreesen, who wrote the first web browser Mosaic that led to the explosive growth of the internet. He has plenty of insightful things to say.  I remember first seeing Mosaic in 1994 as a postdoc in Boulder, Colorado. There I was, doing research that involved programming in C and C++. I was not really happy with what I was doing. I was having a hard time finding the next job. I was one of the first to play around with HTML, and it never occurred to me once that I could pack my bags, move to Silicon Valley, and try to get involved in the burgeoning tech revolution. It just makes me wonder what other obvious things I’m missing right now.

Addendum, 2014-6-5:  Actually, it may have been 1993 that I first saw Mosaic.

## Integrated Information Theory

June 2, 2014

Neuroscientist Giulio Tononi has proposed that consciousness is integrated information and can be measured by a quantity called $\phi$, which is a measure of the amount of information that involves the entire system as a whole. I have never really found this theory to be entirely compelling. While I think that consciousness probably does require some amount of integrated information, I am skeptical that it is the only relevant measure. See here and here for some of my previous thoughts on the topic. One of the reasons that Tononi has proposed a single measure is because it is a way to sidestep what is known as “the hard problem of consciousness”. Instead of trying to explain how a collection of neurons would be endowed with a sense of self-awareness, he posits that consciousness is a property of information and the more $\phi$ one has, the more conscious you become. So in this theory, rocks are not conscious but thermostats are minimally conscious.

Theoretical computer scientist Scott Aaronson has now weighed in on the topic (see here and here). In his inimitable style, Aaronson shows essentially that a large grid of XOR gates could have arbitrarily large $\phi$ and hence be even more conscious than you or me.  He finds this to be highly implausible. Tononi then produced a 14 page response where he essentially doubles down on IIT and claims that indeed a planar array of XOR gates is conscious and we should not be surprised it is so. Aaronson also proposes that we try to solve the “pretty hard problem of consciousness”, which is to come up with a theory or means for deciding when something has consciousness. To me, the fact that we can’t come up with an empirical way to tell whether something is conscious is the best argument for dualism we have. It may even be plausible that the PHPC is undecidable in that solving it would entail the solution of the halting problem. I agree with philosopher David Chalmers (see here) that there are only two possible consistent theories of consciousness. The first is that it is an emergent property of the brain but it has no “causal influence” on events. In other words, consciousness is an epiphenomenon that just allows “us” to be an audience for the dynamical evolution of the universe. The second is that we live in a dualistic world of mind and matter. It is definitely worth reading the posts and the comments, where Chalmers chimes in.

## Did microbes cause the Great Dying?

May 24, 2014

In one of my very first posts almost a decade ago, I wrote about the end-Permian extinction 250 million years ago, which was the greatest mass extinction thus far. In that post I covered research that had ruled out an asteroid impact and found evidence of global warming, possibly due to volcanos, as a cause. Now, a recent paper in PNAS proposes that a horizontal gene transfer event from bacteria to archaea may have been the main cause for the increase of methane and CO2. This paper is one of the best papers I have read in a long time, combining geological field work, mathematical modeling, biochemistry, metabolism, and evolutionary phylogenetic analysis to make a compelling argument for their hypothesis.

Their case hinges on several pieces of evidence. The first comes from well-dated carbon isotopic records from China.  The data shows a steep plunge in the isotopic ratio (i.e ratio between the less abundant but heavier carbon 13 and the lighter more abundant carbon 12) in the inorganic carbonate reservoir with a moderate increase in the organic reservoir. In the earth’s carbon cycle, the organic reservoir comes from the conversion of atmospheric CO2 into carbohydrates via photosynthesis, which prefers carbon 12 to carbon 13. Organic carbon is returned to inorganic form through oxidation by animals eating photosynthetic organisms or by the burning of stored carbon like trees or coal. A steep drop in the isotopic ratio means that there was an extra surge of carbon 12 into the inorganic reservoir. Using a mathematical model, the authors show that in order to explain the steep drop, the inorganic reservoir must have grown superexponentially (faster than exponential). This requires some runaway positive feedback loop that is difficult to explain by geological processes such as volcanic activity, but is something that life is really good at.

The increased methane would have been oxidized to CO2 by other microbes, which would have lowered the oxygen concentration. This would allow for more efficient fermentation and thus more acetate fuel for the archaea to make more methane. The authors showed in another simple mathematical model how this positive feedback loop could lead to superexponential growth. Methane and CO2 are both greenhouse gases and their increase would have caused significant global warming. Anaerobic methane oxidation could also lead to the release of poisonous hydrogen sulfide.

They then considered what microbe could have been responsible. They realized that during the late Permian, a lot of organic material was being deposited in the sediment. The organic reservoir (i.e. fossil fuels, methane hydrates, soil organic matter, peat, etc) was much larger back then than today, as if someone or something used it up at some point. One of the end products of fermentation of this matter would be acetate and that is something archaea like to eat and convert to methane. There are two types of archaea that can do this and one is much more efficient than the other at high acetate concentrations. This increased efficiency was also shown recently to have arisen by a horizontal gene transfer event from a bacterium. A phylogenetic analysis of all known archaea showed that the progenitor of the efficient methanogenic one likely arose 250 million years ago.

The final piece of evidence is that the archaea need nickel to make methane. The authors then looked at the nickel concentrations in their Chinese geological samples and found a sharp increase in nickel immediately before the steep drop in the isotopic ratio. They postulate that the source of the nickel was the massive Siberian volcano eruptions at that time (and previously proposed as the cause of the increased methane and CO2).

This scenario required the unlikely coincidence of several events –  lots of excess organic fuel, low oxygen (and sulfate), increased nickel, and a horizontal gene transfer event. If any of these were missing, the Great Dying may not have taken place. However, given that there have been only 5 mass extinctions, although we may be currently inducing the 6th, low probability events may be required for such calamitous events. This paper should also give us some pause about introducing genetically modified organisms into the environment. While most will probably be harmless, you never know when one will be the match that lights the fire.

## What is the difference between math, science and philsophy?

May 16, 2014

I’ve been listening to the Philosophy Bites podcast recently. One from a few years ago consisted of answers from philosopher’s to the question posed on the spot and without time for deep reflection: What is Philosophy? Some managed to give precise answers, but many struggled. I think one source of conflict they faced as they answered was that they didn’t know how to separate the question of what philosophers actually do from they should be doing. However, I think that a clear distinction between science, math and philosophy as methodologies can be specified precisely. I also think that this is important because practitioner’s in each subject should be aware of what methodology they are actually using and what is appropriate for whatever problem they are working on.

Here are my definitions: Math explores the consequences of rules or assumptions, science is the empirical study of measurable things, and philosophy examines things that cannot be resolved by mathematics or empiricism. With these definitions, practitioner’s of any discipline may use either math, science, or philosophy to help answer whatever question they may be addressing. Scientists need mathematics to work out the consequences of their assumptions and philosophy to help delineate phenomena. Mathematicians need science and philosophy to provide assumptions or rules to analyze. Philosophers need mathematics to sort out arguments and science to test hypotheses experimentally.

Those skeptical of philosophy may suggest that anything that cannot be addressed by math or science has no practical value. However, with these definitions, even the most hardened mathematician or scientist may be practicing philosophy without even knowing it. Atheists like Richard Dawkins should realize that part of their position is based on philosophy and not science. The only truly logical position to take with respect to God is agnosticism. It may be probable that there is not a God that intervenes directly in our lives and that probability may be high but it is not a provable fact. To be an atheist is to put some cutoff on the posterior probability for the existence of God and that cutoff is based on philosophy not science.

While most scientists and mathematicians are cognizant that moral issues may be pertinent to their work (e.g. animal experimentation), they may be less cognizant of what I believe is an equally important philosophical issue , which is the ontological question. Ontology is a philosophical term for the study of what exists. To many pragmatically minded people, this may sound like an ethereal topic (or worse adjective) that has no place in the hard sciences. However, as I pointed out in an earlier post, we can put labels on at most a countably infinite number of things out of an uncountable number of possibilities and for most purposes, our ontological list of things is finite. We thus have to choose and although some of these choices are guided by how we as human agents interact with the world, others will be arbitrary. Determining ontology will involve aspects of philosophy, science and math.

Mathematicians face the ontological problem daily when they decide on what areas to work in and what theorems to prove. The possibilities in mathematics are infinite so it is almost certain that if we were to rerun history some if not many fields would not be reinvented. While scientists may have fewer degrees of freedom to choose from they are also making choices and these choices tend to be confined by history. The ontological problem shows up anytime we try to define a phenomenon. The classification of cognitive disorders is a pure exercise in ontology. Authors of the DSM IV have attempted to be as empirical and objective as possible but there is still plenty of philosophy in their designations of psychiatric conditions. While most string theorists accept that their discipline is mostly mathematical, they should also realize that it is very philosophical. A theory of everything includes the ontology by definition.

Subjects traditionally within the realm of philosophy also have mathematical and scientific aspects. Our morals and values have certainly been shaped by evolution and biological constraints. We should completely rethink our legal philosophy based on what we now know about neuroscience (e.g. see here). The same goes for any discussion of consciousness, the mind-body problem, and free will. To me the real problem with free will isn’t whether or not it exists but rather who or what exactly is exercising that free will and this can be looked at empirically.

So next time when you sit down to solve a problem, think about whether it is one of mathematics, science or philosophy.