Physics of Harrison’s H1 clock

Today’s post is about a subject I presented a few years back at the entrance exams of the french engineering schools. If you like clocks, boats, adventures and unnecessary complications, this is for you 🙂

Historical notice

During the Renaissance the English realised how good they were on seas, and how much potential there was for some ass-kicking with all these loaded European ships making round trips to America. A very crucial problem in sailing is to know where you are. For that you need

  1. Your latitude : how far north are you from the equator ?
  2. Your longitude : how far west are you from London ?

Getting the latitude is easy:

  • During the day, look at the sun (especially at noon). The higher it is, the nearer you are from the equator
  • At night, look at a pole star : the higher it is, the farther away you are from the equator.

But no such tricks to estimate the longitude existed in the 1700’s. Astronomers had thought of something but it was all very complicated and stuff. So the English did what you do when you are out of ideas: they made a law (the Longitude Act) proposing a huge prize (that they had no intention of giving up too easily) to any subject of his Majesty who would come up first with a method.
That’s were John Harrison comes in, with a beautifully simple idea:

  1. Take a clock indicating London’s time onboard.
  2. During your travel, look at the sun to get the time of the place where you are, and compare it to the clock’s time.
  3. If for instance the sun indicates 12h00, and the clock indicates that it is 10h00 in London, then you know that you are \frac{12-10}{24}*360 = 30 degrees west from London.

The only problem is that, if you take your good old pendulum clock on a boat, it will instantly go mad due to all that pitching, rolling and yawing ! So you need to design a more sophisticated device. But John Harrison had a slight advantage here: he was a professional clockmaker !

John Harrison (1693-1776) Credit: Wikipedia

Shake it baby ! Harrison’s H1 clock

In 1730, Harrison imagines his first marine clock, that he calls H1:

Harrison’s H1. Credit : National Maritime Museum, (Greenwich)

Here is an idealized scheme of the beast’s arms :

I believe that the contact between the two arms is there to make sure that they will always move in symetrical ways (otherwise the movement could get chaotic). A very nice feature of this design is that any force that is applied equally to each of the weights is cancelled out and has no effect on the oscillations !

This makes the clock’s movement independent from gravity (meaning that the clock doesn’t need to be vertical to oscillate properly) and from any translational movement of the boat !
However (and that’s where my study started) I noticed that the centrifuge force induced by the rolling of the ship is not compensated, as the weights are not equally distant from the rotation center:

And my (very fundamental and of much historical importance) question was : how much does the rolling affect the clock ?

A small study

Before getting to the real question, let us study how this clock works when it is not perturbated.
Suppose that the left arm, of radius R makes an angle \theta with its equilibrium position (that we suppose to be vertical):

The other arm is in a symetrical position, so the upper spring is streched by a length 2R\sin(\theta) and thus pulls the upper weight of the left arm with a force equal to 2kR\sin(\theta) (where k is the spring’s rigidity constant). The lower spring is compressed by a length 2R\sin(\theta) and the lower weight undergoes a push of 2kR\sin(\theta). These two forces apply negative torques on the arm, both with a lever equal to R\cos(\theta). Finally, neglecting the mass of the arm’s sticks, the arm’s angular momentum is 2R^2M\dot{\theta} where M is the mass of one weight. Once you have all that, the angular momentum theorem gives

2R^2M\ddot{\theta} = -2Rk\sin(\theta)R\cos(\theta)

which after simplifications gives

\ddot{\theta} + \dfrac{k}{M}\sin(2\theta) = 0

Note that (partly because my model is very simple) the oscillations do not depend on the radius of the arms. If the amplitude of the oscillations if really small, we have \sin(2\theta) \simeq 2\theta and the equation becomes

\ddot{\theta} + 2\dfrac{k}{M}\theta = 0

where we recognize an harmonic oscillator of period T = \pi \sqrt{\dfrac{2M}{k}} . If the oscillations are bigger, one can still simulate the dynamics (here with Python):

from pylab import *
from scipy.integrate import odeint

def H1diff(Y,t,k,M):
    Equations describing the dynamics of the model
    Y is [A,dA/dt], t is the current time, k is the spring's constant,
    M is the mass of the weights.
    Returns [dA/dt, ddA/ddt] according to
    dd(A)/ddt + k/M sin(2A) = 0

    A,dA = Y
    return array([dA, -k / M * sin(2*A)])

k=1.0 ; M=1.0 ; initialAngle = pi/6 # parameters
Y0= [pi/6, 0.0] # initial angle pi/6, initial speed 0
ts=linspace(0,10, 100) # 100 times points between 0 and 10
angles = odeint(H1diff,[initialAngle,0],ts,args=(k,M))[:,0] # solving

# plotting
fig,ax = subplots(1)
ax.plot(ts,angles) # plotting
ax.set_xlabel('time',fontsize = 24)
ax.set_ylabel('angle',fontsize = 24)

Or better 🙂

def makeAnimation(k=1.0,M=1.0,nframes=40,R=1.0,**kwargs):
    """ simulates and creates a GIF of the H1 """

    # default for tmax is one period.
    tmax = kwargs.pop('tmax',2*pi*sqrt(M/k))

    # delay is computed so that the animation actually lasts tmax.
    delay = kwargs.pop('delay',int(100.0*tmax/nframes))

    O1 = array([-1.0,0]) # center of the left pendulum
    O2 = array([1.0,0])  # center of the right pendulum
    R = 1.0

    def drawH1(angle,ax):
        """ draws an Harrison H1 clock given the angle of the arms """

        # compute the coordinates of the four weights:

        s,c = sin(angle), cos(angle)
        LU = O1 + R*array([-s,c])
        LD = O1 + R*array([s,-c])
        RU = O2 + R*array([s,c])
        RD = O2 + R*array([-s,-c])

        # draw the springs

        for w1,w2 in (LU,RU),(LD,RD):
            ax.add_line(Line2D(*zip(w1,w2), linewidth=2,color='b'))

        # draw the sticks

        for w1,w2 in (LU,LD),(RU,RD):
            ax.add_line(Line2D(*zip(w1,w2), linewidth=4,color='k'))

        # draw the joints

        for point in O1,O2:
            ax.add_artist(Circle(point,radius = 0.1, color='grey'))

        # draw the weights

        for point in LU,LD,RU,RD:
            ax.add_artist(Circle(point,radius = 0.2,color='k'))

    Y0 = [pi/6,0]
    angles = odeint(H1diff,Y0,ts,args=(k,M))[:,0].flatten()

    fig,ax = subplots(1)

    for i,angle in enumerate(angles):

    # make the gif

    os.system('convert -delay %f -loop -1 hhh*.jpeg harrison.gif'%delay)

    # remove the temporary files
    for myFile in os.listdir('./'):
        if myFile.startswith('hhh'):

H1 simulation

Any role for the Rolling ?

Back to our previous question: how much does the rolling affect the clock’s movement ? Let us just put a few things straight first:

  • It’s impossible to come even close to a solution that I could defend more than one minute.
  • Anyway, who cares ! It’s not like if it was an important problem.

This being said, here is how I went with the modelling. For the clock:

  1. I emailed the Greenwich Museum, where the original H1 is, and they were nice enough to give me a few measures of the clock (length of the arms, radius and material of the weights)
  2. I estimated the maximal angle and checked the period using a video found on the internet.
  3. I then determined the constant k of the springs by simulation and optimisation, as the optimal constant to get a period of one second.

Here is how you do the last step with Python:

from scipy.optimize import fmin

def find_k(T,M,initialAngle):
    """ finds the optimal spring rigidity constant k to obtain an H1
        clock with a period of T given the mass M of the weights
        and the initial angle of the clock """

    def f(k):
        res = odeint(H1diff,[initialAngle,0],ts,args=(k,M))
        finalAngle = res[-1,0]
        return (finalAngle-initialAngle)**2

    k_init = 2*(pi**2)*M / T**2  # expected for very small oscillations
    return fmin( f ,k_init)

To model the boat:

    1. I found a plan on the internet of the very famous HMS Bounty, which is approximately the kind of ship used in Harrison’s times
  1. Since the plan was accurate enough to make a mesh out of it, I used pens and rulers to get the coordinates of the points of the slices of the ship on the plan. then I used some sort of finite elements method (I had not heard of this at the time) to simulate the oscillations of the boat on still water. This gave me a profile of the ship’s rolling over time.
  2. A cool thing : I found on the internet that some people in Australia had made an exact replicate of the Bounty. So I emailed the captain to see if he could measure how long the ship takes to do one roll. He found exactly the same period, more or less one second !

Putting everything together I was able to see how much the rolling impacts the period of the clock over a few minutes, then how much difference it makes in the estimated time after a few days, and how much this represents in terms of distance… Results that I don’t really remember, but once again, who cares ? Isn’t the most important having fun, squeezing equations and poking people from all around the world ?

Epilogue (with spoilers)

The H1 wasn’t precise enough to compute the longitude as precisely as demanded by the Longitude Act. Further improvements kept Harrison busy for his whole life. He produced three clocks, H1, H2, and H3, until he realised that he had been thinking too big, and that the solution might come from a bleeding edge technology of the time: the chronometer ! With its small size and it’s reliable spiral spring, this ancester of the pocket watch was accurate enough to finally get Harrison his prize, after more than 40 years of research !

Finding a subnetwork with a given topology in E. Coli

A few weeks ago I posted about how easy it was, using Python’s Networkx module and the RegulonDB database, to get and navigate through the genetic network of E. Coli. I didn’t mention at the time that you can also find all the subnetworks of E. Coli sharing a given topology with just a few lines of code:

import networkx.algorithms.isomorphism as iso
def find_pattern(graph,pattern,sign_sensitive = True):

    if sign_sensitive:
        edge_match = lambda e1,e2 : (e1['sign']==e2['sign'])
        edge_match = None
    matcher = iso.DiGraphMatcher(graph,pattern,edge_match=edge_match)
    return list(matcher.subgraph_isomorphisms())

Feedforward loops in E. coli

As an example, let us have a look at the feedforward loops in E. coli. A feedforward loop can be described as a gene having an action on another gene, both directly and through the intermediary of another gene, like in the following sketch, where the arrows can represent activations or repressions

A feedforward !

If you want to print all the feedforward loops in E. coli’s network (as represented by RegulonDB), try this :

import networkx as nx
feedforward = nx.DiGraph()
ffwd_in_ecoli = find_pattern(ECN,feedforward,sign_sensitive=False)
for subgraph in ffwd_in_ecoli:
    print subgraph

This prints each feedforward loop with the respective roles of the different genes (A,I, or C). As you can see there are 63 feedforward loops in regulonDB’s network, which is many. Let us now get sign-specific :

import itertools as itt
labels = []
frequencies = []
for signs in itt.product(('+','-'),repeat = 3):

    ffwd = nx.DiGraph()
    frequency = len(list(find_pattern(ECN,ffwd,True)))
    labels.append(" ".join(signs))

frequencies = array([float(f) for f in frequencies])/sum(frequencies)
fig,ax = subplots(1, figsize=(5,5))
title('feedforward loops in E. coli, \n by signs of A->I,I->B,A->B')
pie(frequencies, labels=labels, autopct='%1.1f%%', shadow=True)

Relative frequencies of the different feedforward loops in E. coli

I’m not the first one to do such reasearch, but that won’t prevent me from commenting 🙂 The most frequent pattern can be seen as a double repression : gene A represses gene B, and, as an additional punishment, also represses the activator I of B. Such feedforwards are called coherent, as the two actions of the gene A have the same effect on the gene B. The second most frequent feedforward is incoherent: gene A activates gene B while activating a repressor of B, which seems a little foolish, but can be a way of creating bumps of activity (B is awaken just a few minutes before its repressor puts it down again).

A word of caution

Note that the algorithm didn’t find any feedforward loop with activations only (+++), while some have been reported in the literature. The same way, I couldn’t find any mutual interaction between two genes (A has an action on B and B has an action on A). So if you are going to use regulonDB’s network of interactions for reasearch purposes, always have in mind :

  • The database is INCOMPLETE, meaning that you can use it (carefully) to deduce the existence of stuff, but not to prove the absence of something. In particular, many gene interactions that occur through the intermediary of metabolites are omited ! For instance, the gene cyaA produces the cyclic AMP that is necessary to activate the genes regulated by crp. However only crp appears as a regulator of these genes.
  • In the context of synthetic biology, the few genes added to the bacteria are often supposed to work independently of the rest of the bacteria (some biologists say that they are orthogonal to the rest of the system). Thus, if you design a feedforward loop, it will be a feedforward loop and nothing else. But keep in mind that when the pattern-matching algorithm finds a feedforward loop in E. coli’s networks, each gene of the triangle could also be under the influence of many other, and maybe the behavior of the circuit is not what could be inferred by looking only at these three genes.
  • The algorithm consider isomorphism between subnetworks, which implies that it won’t mind if the subnetwork is under the influence of external genes, but it will mind if the subnetwork has auto-regulations that do not appear on the provided pattern. For instance, fis regulates crp and vice-versa, but this won’t be reported if you simply look for the pattern A\leftrightarrow B because, in addition, fis and crp regulate themselves !

Overstressed ? Here, take some bacteria !

The European jamboree of the 2012 iGEM competition (a synthetic biology contest for undergraduates) is over. And our team didn’t make it through to the world finals at Boston’s MIT 😦 .
Given that the event took place in Amsterdam, maybe this small project of mine would have given us more chances:

Hallucicoli, symbiosis meets happiness

Not sure, however, that it would be appreciated in Boston, as the FBI is one of iGEM’s main sponsors !

So what do you think… Feasible ? Marketable ?

It is a little akward to put such things on the internet, expecially if you, reader, are my future employer, so here is a disclaimer: this is just for fun (note, however, that right now real people are leading a real project to make bacteria produce THC). I do not support the use of drugs under any form and I never, ever had weed in my life (which may not be the case of many of the European iGEMers, who spent three very happy days in Amsterdam !)