text:
# Fun with Rails - Bookmark Manager

I made some time to sit down and play with Rails this evening, and thought I'd work on a simple Bookmark Manager (ala del.icio.us) - a piece of software that I am sorta familiar with, since I use del.icio.us a fair bit, and have some ideas that I think would improve the service. So, I'll try mocking them up in Rails to give me something to chew on while I learn...

So, about half an hour into it, and I have a barebones app that lets me add new bookmarks (complete with title, description, rating, date added, url, etc...), and list existing bookmarks.

Next, to try to work in the concept of Users, and tie bookmarks to specific users. Then, to try to do something interesting with the connections that can happen when multiple users bookmark the same resource, or similar resources, or tag them differently, or give them different ratings, etc...

I'm really liking the concept of Scaffolding in Rails. It feels very similar to starting a WebObjects application via DirectToWeb, and slowly overriding the default functionality as you need to. I also really like the strong MVC separation. Controller files (.rb - analogous to the WebObjects .java files) contain the business logic, while View files (.rhtml - analogous to the WebObjects .html and .wod files merged?) contain HTML with bindings to the methods and variables in the controller.

Almost forgot - I only "wrote" about 3 lines of code. Still, there were several files automagically generated by Rails, and a couple of the .rhtml files with html and bindings, but still - very little actual code to write/manage.

I think I'm really going to like Ruby, as well. It's easily readable without being annoyingly verbose. It's nicely object oriented, without being too abstract. It's interpreted, so it's easy to see code "live" as you tweak it.

I'll play with this app over the next few days/weeks to see what I can come up with. If it gets stable enough, I might deploy it on my desktop box back in the office...

Here's the basic database schema I'm working toward. I still have to think out the bookmark --> tag relationship. I want Tags to be primitive "nouns" - not belonging to any user, and being linked to from Bookmarks through a join table. Adding tags is just a matter of ensure the "noun" exists, and creating (or removing) a join to it.

Update: Thanks to a tip from Hunter, I checked out the acts_as_taggable plugin/extension/whatever for Rails (description here). Basically, you just drop this bad boy into place, and declare in your controller that the class acts_as_taggable - then, the plugin takes care of managing a many-to-many relationship between that class and the tags.

Sample code from the description page for acts_as_taggable:

class Photo < ActiveRecord::Base
acts_as_taggable
end

elephant = Photo.find(4437)
elephant.tag 'zoo animals nature'

elephant.tagged_with?('urban') # => false
elephant.tags.size # => 3
elephant.tag_names # => [ 'zoo', 'animals', 'nature' ]

Can it get any easier than that? It's doing exactly what I was describing as wanting, but For Free. Sweet.

Update 2: Yes, it can get easier. And more powerful. The developer of acts_as_taggable reworked the mixin/plugin after some thought, and came up with something even cooler (documentation here).

In poking around some of the Rails sites, I keep thinking that although it's still a bit immature (it just hit 1.0 a month ago), it "feels" like WebObjects, but based on Open Source from soup to nuts, and with an apparently vibrant and cool developer community. WebObjects is waaaaay more mature, but isn't Open Source (hence the licensing fiasco this summer), and it feels like the developer community may be atrophying a bit (although I could be wrong on that - I hope I am).

Update 3: About an hour into noodling with the acts_as_taggable mixin, and I have a fully folksonomic bookmark manager. I've got a simple form that lets me add bookmarks, complete with a freeform text field for entering tags, which get bursted into individual Tag rows with the appropriate join table invisibly managed.

I've got some more work to do to clean it up, and then I want to add in the concept of users - which apparently is absolutely trivial in Rails. Then, add in some OPML and RSS exports, and support for bookmarklets, and perhaps implementing part of the del.icio.us API on the back end to allow other things to talk to it, and I've got a decent bookmark manager, with only a handful of code written. Sweet.


text:
# Spiking Neural Systems

This is an introduction to spiking neural systems with Julia's DifferentialEquations package. We will cover three different models: leaky integrate-and-fire, Izhikevich, and Hodgkin-Huxley. Finally we will also learn about two mechanisms that simulate synaptic inputs like real neurons receive them. The alpha synapse and the Tsodyks-Markram synapse. Let's get started with the leaky integrate-and-fire (LIF) model.

## The Leaky Integrate-and-Fire Model

The LIF model is an extension of the integrate-and-fire (IF) model. While the IF model simply integrates input until it fires, the LIF model integrates input but also decays towards an equilibrium potential. This means that inputs that arrive in quick succession have a much higher chance to make the cell spike as opposed to inputs that are further apart in time. The LIF is a more realistic neuron model than the IF, because it is known from real neurons that the timing of inputs is extremely relevant for their spiking.

The LIF model has five parameters, gL, EL, C, Vth, I and we define it in the lif(u, p, t) function.

using DifferentialEquations
using Plots
gr()

function lif(u,p,t);
gL, EL, C, Vth, I = p
(-gL*(u-EL)+I)/C
end
lif (generic function with 1 method)

Our system is described by one differential equation: (-gL*(u-EL)+I)/C, where u is the voltage, I is the input, gL is the leak conductance, EL is the equilibrium potential of the leak conductance and C is the membrane capacitance. Generally, any change of the voltage is slowed down (filtered) by the membrane capacitance. That's why we divide the whole equation by C. Without any external input, the voltage always converges towards EL. If u is larger than EL, u decreases until it is at EL. If u is smaller than EL, u increases until it is at EL. The only other thing that can change the voltage is the external input I.

Our lif function requires a certain parameter structure because it will need to be compatible with the DifferentialEquations interface. The input signature is lif(u, p, t) where u is the voltage, p is the collection of the parameters that describe the equation and t is time. You might wonder why time does not show up in our equation, although we need to calculate the change in voltage with respect to time. The ODE solver will take care of time for us. One of the advantages of the ODE solver as opposed to calculating the change of u in a for loop is that many ODE solver algorithms can dynamically adjust the time step in a way that is efficient and accurate.

One crucial thing is still missing however. This is supposed to be a model of neural spiking, right? So we need a mechanism that recognizes the spike and hyperpolarizes u in response. For this purpose we will use callbacks. They can make discontinuous changes to the model when certain conditions are met.

function thr(u,t,integrator)
integrator.u > integrator.p[4]
end

function reset!(integrator)
integrator.u = integrator.p[2]
end

threshold = DiscreteCallback(thr,reset!)
current_step= PresetTimeCallback([2,15],integrator -> integrator.p[5] += 210.0)
cb = CallbackSet(current_step,threshold)
DiffEqBase.CallbackSet{Tuple{}, Tuple{DiffEqBase.DiscreteCallback{DiffEqCal
lbacks.var"#61#64"{Vector{Int64}}, DiffEqCallbacks.var"#62#65"{Main.##Weave
SandBox#10480.var"#1#2"}, DiffEqCallbacks.var"#63#66"{typeof(DiffEqBase.INI
TIALIZE_DEFAULT), Bool, Vector{Int64}, Main.##WeaveSandBox#10480.var"#1#2"}
, typeof(DiffEqBase.FINALIZE_DEFAULT)}, DiffEqBase.DiscreteCallback{typeof(
Main.##WeaveSandBox#10480.thr), typeof(Main.##WeaveSandBox#10480.reset!), t
ypeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}
}((), (DiffEqBase.DiscreteCallback{DiffEqCallbacks.var"#61#64"{Vector{Int64
}}, DiffEqCallbacks.var"#62#65"{Main.##WeaveSandBox#10480.var"#1#2"}, DiffE
qCallbacks.var"#63#66"{typeof(DiffEqBase.INITIALIZE_DEFAULT), Bool, Vector{
Int64}, Main.##WeaveSandBox#10480.var"#1#2"}, typeof(DiffEqBase.FINALIZE_DE
FAULT)}(DiffEqCallbacks.var"#61#64"{Vector{Int64}}([2, 15]), DiffEqCallback
s.var"#62#65"{Main.##WeaveSandBox#10480.var"#1#2"}(Main.##WeaveSandBox#1048
0.var"#1#2"()), DiffEqCallbacks.var"#63#66"{typeof(DiffEqBase.INITIALIZE_DE
FAULT), Bool, Vector{Int64}, Main.##WeaveSandBox#10480.var"#1#2"}(DiffEqBas
e.INITIALIZE_DEFAULT, true, [2, 15], Main.##WeaveSandBox#10480.var"#1#2"())
, DiffEqBase.FINALIZE_DEFAULT, Bool[1, 1]), DiffEqBase.DiscreteCallback{typ
eof(Main.##WeaveSandBox#10480.thr), typeof(Main.##WeaveSandBox#10480.reset!
), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAUL
T)}(Main.##WeaveSandBox#10480.thr, Main.##WeaveSandBox#10480.reset!, DiffEq
Base.INITIALIZE_DEFAULT, DiffEqBase.FINALIZE_DEFAULT, Bool[1, 1])))

Our condition is thr(u,t,integrator) and the condition kicks in when integrator.u > integrator.p[4] where p[4] is our threshold parameter Vth. Our effect of the condition is reset!(integrator). It sets u back to the equilibrium potential p[2]. We then wrap both the condition and the effect into a DiscreteCallback called threshold. There is one more callback called PresetTimeCallback that is particularly useful. This one increases the input p[5] at t=2 and t=15 by 210.0. Both callbacks are then combined into a CallbackSet. We are almost done to simulate our system we just need to put numbers on our initial voltage and parameters.

u0 = -75
tspan = (0.0, 40.0)
# p = (gL, EL, C, Vth, I)
p = [10.0, -75.0, 5.0, -55.0, 0]

prob = ODEProblem(lif, u0, tspan, p, callback=cb)
ODEProblem with uType Int64 and tType Float64. In-place: false
timespan: (0.0, 40.0)
u0: -75

Our initial voltage is u0 = - 75, which will be the same as our equilibrium potential, so we start at a stable point. Then we define the timespan we want to simulate. The time scale of the LIF as it is defined conforms roughly to milliseconds. Then we define our parameters as p = [10.0, -75.0, 5.0, -55.0, 0]. Remember that gL, EL, C, Vth, I = p. Finally we wrap everything into a call to ODEProblem. Can't forget the CallbackSet. With that our model is defined. Now we just need to solve it with a quick call to solve.

sol = solve(prob)
retcode: Success
Interpolation: automatic order switching interpolation
t: 153-element Vector{Float64}:
0.0
9.999999999999999e-5
0.0010999999999999998
0.011099999999999997
0.11109999999999996
1.1110999999999995
2.0
2.0
2.6300346673750097
2.9226049547524595
⋮
38.34157935968204
38.78215179003683
38.78215179003683
39.222724173706894
39.222724173706894
39.6632965982261
39.6632965982261
40.0
40.0
u: 153-element Vector{Float64}:
-75.0
-75.0
-75.0
-75.0
-75.0
-75.0
-75.0
-75.0
-59.978080111690375
-57.32999167299642
⋮
-75.0
-50.40489310815222
-75.0
-50.404894730067554
-75.0
-50.404893310891545
-75.0
-54.419318668318546
-75.0

First of all the solve output tells us if solving the system generally worked. In this case we know it worked because the return code (retcode) says Success. Then we get the numbers for the timestep and the solution to u. The raw numbers are not super interesting to let's plot our solution.

plot(sol)

We see that the model is resting at -75 while there is no input. At t=2 the input increases by 210 and the model starts to spike. Spiking does not start immediately because the input first has to charge the membrane capacitance. Notice how once spiking starts it very quickly becomes extremely regular. Increasing the input again at t=15 increases firing as we would expect but it is still extremely regular. This is one of the features of the LIF. The firing frequency is regular for constant input and a linear function of the input strength. There are ways to make LIF models less regular. For example we could use certain noise types at the input. We could also simulate a large number of LIF models and connect them synaptically. Instead of going into those topics, we will move on to the Izhikevich model, which is known for its ability to generate a large variety of spiking dynamics during constant inputs.

## The Izhikevich Model

The Izhikevich model is a two-dimensional model of neuronal spiking. It was derived from a bifurcation analysis of a cortical neuron. Because it is two-dimensional it can generate much more complex spike dynamics than the LIF model. The kind of dynamics depend on the four parameters and the input a, b, c, d, I = p. All the concepts are the same as above, expect for some minor changes to our function definitions to accomodate for the second dimension.

#Izhikevichch Model
using DifferentialEquations
using Plots

function izh!(du,u,p,t);
a, b, c, d, I = p

du[1] = 0.04*u[1]^2+5*u[1]+140-u[2]+I
du[2] = a*(b*u[1]-u[2])
end
izh! (generic function with 1 method)

This is our Izhikevich model. There are two important changes here. First of all, note the additional input parameter du. This is a sequence of differences. du[1] corresponds to the voltage (the first dimension of the system) and du[2] corresponds to the second dimension. This second dimension is called u in the original Izhikevich work amnd it makes the notation a little annoying. In this tutorial I will generally stick to Julia and DifferentialEquations conventions as opposed to conventions of the specific models and du is commonly used. We will never define du ourselves outside of the function but the ODE solver will use it internally. The other change here is the ! after our function name. This signifies that du will be preallocated before integration and then updated in-place, which saves a lot of allocation time. Now we just need our callbacks to take care of spikes and increase the input.

function thr(u,t,integrator)
integrator.u[1] >= 30
end

function reset!(integrator)
integrator.u[1] = integrator.p[3]
integrator.u[2] += integrator.p[4]
end

threshold = DiscreteCallback(thr,reset!)
current_step= PresetTimeCallback(50,integrator -> integrator.p[5] += 10)
cb = CallbackSet(current_step,threshold)
DiffEqBase.CallbackSet{Tuple{}, Tuple{DiffEqBase.DiscreteCallback{DiffEqCal
lbacks.var"#61#64"{Int64}, DiffEqCallbacks.var"#62#65"{Main.##WeaveSandBox#
10480.var"#3#4"}, DiffEqCallbacks.var"#63#66"{typeof(DiffEqBase.INITIALIZE_
DEFAULT), Bool, Int64, Main.##WeaveSandBox#10480.var"#3#4"}, typeof(DiffEqB
ase.FINALIZE_DEFAULT)}, DiffEqBase.DiscreteCallback{typeof(Main.##WeaveSand
Box#10480.thr), typeof(Main.##WeaveSandBox#10480.reset!), typeof(DiffEqBase
.INITIALIZE_DEFAULT), typeof(DiffEqBase.FINALIZE_DEFAULT)}}}((), (DiffEqBas
e.DiscreteCallback{DiffEqCallbacks.var"#61#64"{Int64}, DiffEqCallbacks.var"
#62#65"{Main.##WeaveSandBox#10480.var"#3#4"}, DiffEqCallbacks.var"#63#66"{t
ypeof(DiffEqBase.INITIALIZE_DEFAULT), Bool, Int64, Main.##WeaveSandBox#1048
0.var"#3#4"}, typeof(DiffEqBase.FINALIZE_DEFAULT)}(DiffEqCallbacks.var"#61#
64"{Int64}(50), DiffEqCallbacks.var"#62#65"{Main.##WeaveSandBox#10480.var"#
3#4"}(Main.##WeaveSandBox#10480.var"#3#4"()), DiffEqCallbacks.var"#63#66"{t
ypeof(DiffEqBase.INITIALIZE_DEFAULT), Bool, Int64, Main.##WeaveSandBox#1048
0.var"#3#4"}(DiffEqBase.INITIALIZE_DEFAULT, true, 50, Main.##WeaveSandBox#1
0480.var"#3#4"()), DiffEqBase.FINALIZE_DEFAULT, Bool[1, 1]), DiffEqBase.Dis
creteCallback{typeof(Main.##WeaveSandBox#10480.thr), typeof(Main.##WeaveSan
dBox#10480.reset!), typeof(DiffEqBase.INITIALIZE_DEFAULT), typeof(DiffEqBas
e.FINALIZE_DEFAULT)}(Main.##WeaveSandBox#10480.thr, Main.##WeaveSandBox#104
80.reset!, DiffEqBase.INITIALIZE_DEFAULT, DiffEqBase.FINALIZE_DEFAULT, Bool
[1, 1])))

One key feature of the Izhikevich model is that each spike increases our second dimension u[2] by a preset amount p[4]. Between spikes u[2] decays to a value that depends on p[1] and p[2] and the equilibrium potential p[3]. Otherwise the code is not too different from the LIF model. We will again need to define our parameters and we are ready to simulate.

p = [0.02, 0.2, -50, 2, 0]
u0 = [-65, p[2]*-65]
tspan = (0.0, 300)

prob = ODEProblem(izh!, u0, tspan, p, callback=cb)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 300.0)
u0: 2-element Vector{Float64}:
-65.0
-13.0
sol = solve(prob);
plot(sol, vars=1)

This spiking type is called chattering. It fires with intermittent periods of silence. Note that the input starts at t=50 and remain constant for the duration of the simulation. One of mechanisms that sustains this type of firing is the spike induced hyperpolarization coming from our second dimension, so let's look at this variable.

plot(sol, vars=2)

Our second dimension u[2] increases with every spike. When it becomes too large, the system cannot generate another spike until u[2] has decayed to a value small enough that spiking can resume. This process repeats. In this model, spiking is no longer regular like it was in the LIF. Here we have two frequencies, the frequency during the spiking state and the frequency between spiking states. The LIF model was dominated by one single frequency that was a function of the input strength. Let's see if we can generate another spiking type by changing the parameters.

p = [0.02, 0.2, -65, 8, 0]
u0 = [-65, p[2]*-65]
tspan = (0.0, 300)

prob = ODEProblem(izh!, u0, tspan, p, callback=cb)
sol = solve(prob);
plot(sol, vars=1)

This type is called regularly spiking and we created it just by lowering p[3] and increasing p[4]. Note that the type is called regularly spiking but it is not instantaneously regular. The instantenous frequency is higher in the beginning. This is called spike frequency adaptation and is a common property of real neurons. There are many more spike types that can be generated. Check out the original Izhikevich work and create your own favorite neuron!

## Hodgkin-Huxley Model

The Hodgkin-Huxley (HH) model is our first biophysically realistic model. This means that all parameters and mechanisms of the model represent biological mechanisms. Specifically, the HH model simulates the ionic currents that depolarize and hyperpolarize a neuron during an action potential. This makes the HH model four-dimensional. Let's see how it looks.

using DifferentialEquations
using Plots

# Potassium ion-channel rate functions
alpha_n(v) = (0.02 * (v - 25.0)) / (1.0 - exp((-1.0 * (v - 25.0)) / 9.0))
beta_n(v) = (-0.002 * (v - 25.0)) / (1.0 - exp((v - 25.0) / 9.0))

# Sodium ion-channel rate functions
alpha_m(v) = (0.182*(v + 35.0)) / (1.0 - exp((-1.0 * (v + 35.0)) / 9.0))
beta_m(v) = (-0.124 * (v + 35.0)) / (1.0 - exp((v + 35.0) / 9.0))

alpha_h(v) = 0.25 * exp((-1.0 * (v + 90.0)) / 12.0)
beta_h(v) = (0.25 * exp((v + 62.0) / 6.0)) / exp((v + 90.0) / 12.0)

function HH!(du,u,p,t);
gK, gNa, gL, EK, ENa, EL, C, I = p
v, n, m, h = u

du[1] = (-(gK * (n^4.0) * (v - EK)) - (gNa * (m ^ 3.0) * h * (v - ENa)) - (gL * (v - EL)) + I) / C
du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)
end
HH! (generic function with 1 method)

We have three different types of ionic conductances. Potassium, sodium and the leak. The potassium and sodium conducance are voltage gated. They increase or decrease depending on the voltage. In ion channel terms, open channels can transition to the closed state and closed channels can transition to the open state. It's probably easiest to start with the potassium current described by gK * (n^4.0) * (EK - v). Here gK is the total possible conductance that we could reach if all potassium channels were open. If all channels were open, n would equal 1 which is usually not the case. The transition from open state to closed state is modeled in alpha_n(v) while the transition from closed to open is in beta_n(v). Because potassium conductance is voltage gated, these transitions depend on v. The numbers in alpha_n; beta_n were calculated by Hodgkin and Huxley based on their extensive experiments on the squid giant axon. They also determined, that n needs to be taken to the power of 4 to correctly model the amount of open channels.

The sodium current is not very different but it has two gating variables, m, h instead of one. The leak conductance gL has no gating variables because it is not voltage gated. Let's move on to the parameters. If you want all the details on the HH model you can find a great description here.

current_step= PresetTimeCallback(100,integrator -> integrator.p[8] += 1)

# n, m & h steady-states
n_inf(v) = alpha_n(v) / (alpha_n(v) + beta_n(v))
m_inf(v) = alpha_m(v) / (alpha_m(v) + beta_m(v))
h_inf(v) = alpha_h(v) / (alpha_h(v) + beta_h(v))

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60)]
tspan = (0.0, 1000)

prob = ODEProblem(HH!, u0, tspan, p, callback=current_step)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 1000.0)
u0: 4-element Vector{Float64}:
-60.0
0.0007906538330645917
0.08362733690208038
0.41742979353768533

For the HH model we need only one callback. The PresetTimeCallback that starts our input current. We don't need to reset the voltage when it reaches threshold because the HH model has its own repolarization mechanism. That is the potassium current, which activates at large voltages and makes the voltage more negative. The three functions n_inf; m_inf; h_inf help us to find good initial values for the gating variables. Those functions tell us that the steady-state gating values should be for the initial voltage. The parameters were chosen in a way that the properties of the model roughly resemble that of a cortical pyramidal cell instead of the giant axon Hodgkin and Huxley were originally working on.

sol = solve(prob);
plot(sol, vars=1)

That's some good regular voltage spiking. One of the cool things about a biophysically realistic model is that the gating variables tell us something about the mechanisms behind the action potential. You might have seen something like the following plot in a biology textbook.

plot(sol, vars=[2,3,4], tspan=(105.0,130.0))

So far we have only given our neurons very simple step inputs by simply changing the number I. Actual neurons recieve their inputs mostly from chemical synapses. They produce conductance changes with very complex structures. In the next chapter we will try to incorporate a synapse into our HH model.

## Alpha Synapse

One of the most simple synaptic mechanisms used in computational neuroscience is the alpha synapse. When this mechanism is triggered, it causes an instantanouse rise in conductance followed by an exponential decay. Let's incorporate that into our HH model.

function gSyn(max_gsyn, tau, tf, t);
if t-tf >= 0
return max_gsyn * exp(-(t-tf)/tau)
else
return 0.0
end
end
function HH!(du,u,p,t);
gK, gNa, gL, EK, ENa, EL, C, I, max_gSyn, ESyn, tau, tf = p
v, n, m, h = u

ISyn = gSyn(max_gSyn, tau, tf, t) * (v - ESyn)

du[1] = (-(gK * (n^4.0) * (v - EK)) - (gNa * (m ^ 3.0) * h * (v - ENa)) - (gL * (v - EL)) + I - ISyn) / C
du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)
end
HH! (generic function with 1 method)

gSyn models the step to the maximum conductance and the following exponential decay with time constant tau. Of course we only want to integrate the conductance at and after time tf, the onset of the synaptic response. Before tf, gSyn returns zero. To convert the conductance to a current, we multiply by the difference between the current voltage and the synapses equilibrium voltage: ISyn = gSyn(max_gSyn, tau, tf, t) * (v - ESyn). Later we will set the parameter ESyn to 0, making this synapse an excitatory synapse. Excitatory synapses have equilibrium potentials far above the resting potential. Let's see what our synapse does to the voltage of the cell.

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 0.008, 0, 20, 100]
tspan = (0.0, 200)
prob = ODEProblem(HH!, u0, tspan, p)
sol = solve(prob);
plot(sol, vars=1)

What you see here is called an excitatory postsynaptic potential (EPSP). It is the voltage response to a synaptic current. While our synaptic conductance rises instantly, the voltage response rises at a slower time course that is given by the membrane capacitance C. This particular voltage response is not strong enough to evoke spiking, so we say it is subthreshold. To get a suprathreshold response that evokes spiking we simply increase the parameter max_gSyn to increase the maximum conductance.

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 0.01, 0, 20, 100]
tspan = (0.0, 200)
prob = ODEProblem(HH!, u0, tspan, p)
sol = solve(prob);
plot!(sol, vars=1)

This plot shows both the subthreshold EPSP from above as well as the suprathreshold EPSP. Alpha synapses are nice because of their simplicity. Real synapses however, are extremely complex structures. One of the most important features of real synapses is that their maximum conductance is not the same on every event. The number and frequency of synaptic events changes the size of the maximum conductance in a dynamic way. While we usually avoid anatomical and biophysical details of real synapses, there is a widely used phenomenological way to capture those dynamics called the Tsodyks-Markram synapse.

## Tsodyks-Markram Synapse

The Tsodyks-Markram synapse (TMS) is a dynamic system that models the changes of maximum conductance that occur between EPSPs at different frequencies. The single response is similar to the alpha synapse in that it rises instantaneously and decays exponentially. The maximum conductance it reaches depends on the event history. To simulate the TMS we need to incorporate three more dimensions, u, R, gsyn into our system. u decays towards 0, R decays towards 1 and gsyn decays towards 0 as it did with the alpha synapse. The crucial part of the TMS is in epsp!, where we handle the discontinuities when a synaptic event occurs. Instead of just setting gsyn to the maximum conductance gmax, we increment gsyn by a fraction of gmax that depends on the other two dynamic parameters. The frequency dependence comes from the size of the time constants tau_u and tau_R. Enough talk, let's simulate it.

function HH!(du,u,p,t);
gK, gNa, gL, EK, ENa, EL, C, I, tau, tau_u, tau_R, u0, gmax, Esyn  = p
v, n, m, h, u, R, gsyn = u

du[1] = ((gK * (n^4.0) * (EK - v)) + (gNa * (m ^ 3.0) * h * (ENa - v)) + (gL * (EL - v)) + I + gsyn * (Esyn - v)) / C
du[2] = (alpha_n(v) * (1.0 - n)) - (beta_n(v) * n)
du[3] = (alpha_m(v) * (1.0 - m)) - (beta_m(v) * m)
du[4] = (alpha_h(v) * (1.0 - h)) - (beta_h(v) * h)

# Synaptic variables
du[5] = -(u/tau_u)
du[6] = (1-R)/tau_R
du[7] = -(gsyn/tau)
end

function epsp!(integrator);
integrator.u[5] += integrator.p[12] * (1 - integrator.u[5])
integrator.u[7] += integrator.p[13] * integrator.u[5] * integrator.u[6]
integrator.u[6] -= integrator.u[5] * integrator.u[6]

end

epsp_ts= PresetTimeCallback(100:100:500, epsp!)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 30, 1000, 50, 0.5, 0.005, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60), 0.0, 1.0, 0.0]
tspan = (0.0, 700)
prob = ODEProblem(HH!, u0, tspan, p, callback=epsp_ts)
sol = solve(prob);
plot(sol, vars=1)

plot(sol, vars=7)

Both the voltage response as well as the conductances show what is called short-term facilitation. An increase in peak conductance over multiple synaptic events. Here the first event has a conductance of around 0.0025 and the last one of 0.004. We can plot the other two varialbes to see what underlies those dynamics

plot(sol, vars=[5,6])

Because of the time courses at play here, this facilitation is frequency dependent. If we increase the period between these events, facilitation does not occur.

epsp_ts= PresetTimeCallback(100:1000:5100, epsp!)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 30, 500, 50, 0.5, 0.005, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60), 0.0, 1.0, 0.0]
tspan = (0.0, 5300)
prob = ODEProblem(HH!, u0, tspan, p, callback=epsp_ts)
sol = solve(prob);
plot(sol, vars=7)

plot(sol, vars=[5,6])

We can also change these time constants such that the dynamics show short-term depression instead of facilitation.

epsp_ts= PresetTimeCallback(100:100:500, epsp!)

p = [35.0, 40.0, 0.3, -77.0, 55.0, -65.0, 1, 0, 30, 100, 1000, 0.5, 0.005, 0]
u0 = [-60, n_inf(-60), m_inf(-60), h_inf(-60), 0.0, 1.0, 0.0]
tspan = (0.0, 700)
prob = ODEProblem(HH!, u0, tspan, p, callback=epsp_ts)
sol = solve(prob);
plot(sol, vars=7)

plot(sol, vars=[5,6])

Just changing those two time constants has changed the dynamics to short-term depression. This is still frequency dependent. Changing these parameters can generate a variety of different short-term dynamics.

## Summary

That's it for now. Thanks for making it this far. If you want to learn more about neuronal dynamics, this is a great resource. If you want to learn more about Julia check out the official website and to learn more about the DifferentialEquations package you are in the right place, because this chapter is part of a larger tutorial series about just that.

## Appendix

These tutorials are a part of the SciMLTutorials.jl repository, found at: https://github.com/SciML/SciMLTutorials.jl. For more information on high-performance scientific machine learning, check out the SciML Open Source Software Organization https://sciml.ai.

To locally run this tutorial, do the following commands:

using SciMLTutorials
SciMLTutorials.weave_file("tutorials/models","08-spiking_neural_systems.jmd")

Computer Information:

Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
OS: Linux (x86_64-pc-linux-gnu)
CPU: AMD EPYC 7502 32-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-11.0.1 (ORCJIT, znver2)
Environment:
JULIA_DEPOT_PATH = /root/.cache/julia-buildkite-plugin/depots/a6029d3a-f78b-41ea-bc97-28aa57c6c6ea

Package Information:

      Status /var/lib/buildkite-agent/builds/1-amdci4-julia-csail-mit-edu/julialang/scimltutorials-dot-jl/tutorials/models/Project.toml
[479239e8] Catalyst v6.12.1
[459566f4] DiffEqCallbacks v2.16.1
[f3b72e0c] DiffEqDevTools v2.27.2
[055956cb] DiffEqPhysics v3.9.0
[0c46a032] DifferentialEquations v6.17.1
[31c24e10] Distributions v0.24.18
[587475ba] Flux v0.12.1
[f6369f11] ForwardDiff v0.10.18
[23fbe1c1] Latexify v0.15.5
[961ee093] ModelingToolkit v5.17.3
[2774e3e8] NLsolve v4.5.1
[315f7962] NeuralPDE v3.10.1
[429524aa] Optim v1.3.0
[1dea7af3] OrdinaryDiffEq v5.56.0
[91a5bcdd] Plots v1.15.2
[731186ca] RecursiveArrayTools v2.11.4
[30cb0354] SciMLTutorials v0.9.0
[789caeaf] StochasticDiffEq v6.34.1
[37e2e46d] LinearAlgebra
[2f01184e] SparseArrays

And the full manifest:

      Status /var/lib/buildkite-agent/builds/1-amdci4-julia-csail-mit-edu/julialang/scimltutorials-dot-jl/tutorials/models/Manifest.toml
[c3fe647b] AbstractAlgebra v0.16.0
[621f4979] AbstractFFTs v1.0.1
[1520ce14] AbstractTrees v0.3.4
[ec485272] ArnoldiMethod v0.1.0
[4fba245c] ArrayInterface v3.1.15
[4c555306] ArrayLayouts v0.7.0
[13072b0f] AxisAlgorithms v1.0.0
[ab4f0b2a] BFloat16s v0.1.0
[aae01518] BandedMatrices v0.16.9
[764a87c0] BoundaryValueDiffEq v2.7.1
[fa961155] CEnum v0.4.1
[00ebfdb7] CSTParser v2.5.0
[052768ef] CUDA v2.6.3
[479239e8] Catalyst v6.12.1
[082447d4] ChainRules v0.7.65
[d360d2e6] ChainRulesCore v0.9.44
[944b1d66] CodecZlib v0.7.0
[35d6a980] ColorSchemes v3.12.1
[3da002f7] ColorTypes v0.11.0
[5ae59095] Colors v0.12.8
[861a8166] Combinatorics v1.0.2
[a80b9123] CommonMark v0.8.1
[38540f10] CommonSolve v0.2.0
[bbf7d656] CommonSubexpressions v0.3.0
[34da2185] Compat v3.30.0
[8f4d0f93] Conda v1.5.2
[88cd18e8] ConsoleProgressMonitor v0.1.2
[187b0558] ConstructionBase v1.2.1
[d38c429a] Contour v0.5.7
[a8cc5b0e] Crayons v4.0.4
[8a292aeb] Cuba v2.2.0
[667455a9] Cubature v1.5.1
[9a962f9c] DataAPI v1.6.0
[82cc6244] DataInterpolations v3.3.1
[864edb3b] DataStructures v0.18.9
[e2d170a0] DataValueInterfaces v1.0.0
[bcd4f6db] DelayDiffEq v5.31.0
[2b5f629d] DiffEqBase v6.62.2
[459566f4] DiffEqCallbacks v2.16.1
[f3b72e0c] DiffEqDevTools v2.27.2
[5a0ffddc] DiffEqFinancial v2.4.0
[aae7a2af] DiffEqFlux v1.37.0
[c894b116] DiffEqJump v6.14.2
[77a26b50] DiffEqNoiseProcess v5.7.3
[055956cb] DiffEqPhysics v3.9.0
[41bf760c] DiffEqSensitivity v6.45.0
[163ba53b] DiffResults v1.0.3
[b552c78f] DiffRules v1.0.2
[0c46a032] DifferentialEquations v6.17.1
[c619ae07] DimensionalPlotRecipes v1.2.0
[b4f34e82] Distances v0.10.3
[31c24e10] Distributions v0.24.18
[ffbed154] DocStringExtensions v0.8.4
[e30172f5] Documenter v0.26.3
[d4d017d3] ExponentialUtilities v1.8.4
[e2ba6199] ExprTools v0.1.3
[c87230d0] FFMPEG v0.4.0
[7a1cc6ca] FFTW v1.4.1
[9aa1b823] FastClosures v0.3.2
[1a297f60] FillArrays v0.11.7
[6a86dc24] FiniteDiff v2.8.0
[53c48c17] FixedPointNumbers v0.8.4
[587475ba] Flux v0.12.1
[59287772] Formatting v0.4.2
[f6369f11] ForwardDiff v0.10.18
[069b7b12] FunctionWrappers v1.1.2
[d9f16b24] Functors v0.2.1
[0c68f7d7] GPUArrays v6.4.1
[61eb1bfa] GPUCompiler v0.10.0
[28b8d3ca] GR v0.57.4
[a75be94c] GalacticOptim v1.2.0
[5c1252a2] GeometryBasics v0.3.12
[af5da776] GlobalSensitivity v1.0.0
[42e2da0e] Grisu v1.0.2
[19dc6840] HCubature v1.5.0
[cd3eb016] HTTP v0.9.9
[eafb193a] Highlights v0.4.5
[0e44f5e4] Hwloc v2.0.0
[7073ff75] IJulia v1.23.2
[b5f81e59] IOCapture v0.1.1
[7869d1d1] IRTools v0.4.2
[615f187c] IfElse v0.1.0
[d25df0c9] Inflate v0.1.2
[83e8ac13] IniFile v0.5.0
[a98d9a8b] Interpolations v0.13.2
[c8e1da08] IterTools v1.3.0
[42fd0dbc] IterativeSolvers v0.9.1
[82899510] IteratorInterfaceExtensions v1.0.0
[692b3bcd] JLLWrappers v1.3.0
[682c06a0] JSON v0.21.1
[98e50ef6] JuliaFormatter v0.13.7
[e5e0dc1b] Juno v0.8.4
[5ab0869b] KernelDensity v0.6.3
[929cbde3] LLVM v3.7.1
[b964fa9f] LaTeXStrings v1.2.1
[2ee39098] LabelledArrays v1.6.1
[23fbe1c1] Latexify v0.15.5
[a5e1c1ea] LatinHypercubeSampling v1.8.0
[73f95e8e] LatticeRules v0.0.1
[093fc24a] LightGraphs v1.3.5
[d3d80556] LineSearches v7.1.1
[2ab3a3ac] LogExpFunctions v0.2.4
[e6f89c97] LoggingExtras v0.4.6
[bdcacae8] LoopVectorization v0.12.23
[1914dd2f] MacroTools v0.5.6
[739be429] MbedTLS v1.0.3
[442fdcdd] Measures v0.3.1
[e89f7d12] Media v0.5.0
[c03570c3] Memoize v0.4.4
[e1d29d7a] Missings v1.0.0
[961ee093] ModelingToolkit v5.17.3
[4886b29c] MonteCarloIntegration v0.0.2
[f9640e96] MultiScaleArrays v1.8.1
[ffc61752] Mustache v1.0.10
[d41bc354] NLSolversBase v7.8.0
[2774e3e8] NLsolve v4.5.1
[872c559c] NNlib v0.7.19
[77ba4419] NaNMath v0.3.5
[315f7962] NeuralPDE v3.10.1
[8913a72c] NonlinearSolve v0.3.8
[6fe1bfb0] OffsetArrays v1.9.0
[429524aa] Optim v1.3.0
[bac558e1] OrderedCollections v1.4.1
[1dea7af3] OrdinaryDiffEq v5.56.0
[90014a1f] PDMats v0.11.0
[65888b18] ParameterizedFunctions v5.10.0
[d96e819e] Parameters v0.12.2
[69de0a69] Parsers v1.1.0
[995b91a9] PlotUtils v1.0.10
[91a5bcdd] Plots v1.15.2
[e409e4f3] PoissonRandom v0.4.0
[f517fe37] Polyester v0.3.1
[85a6dd25] PositiveFactorizations v0.2.4
[21216c6a] Preferences v1.2.2
[33c8b6b6] ProgressLogging v0.1.4
[92933f4c] ProgressMeter v1.6.2
[8a4e6c94] QuasiMonteCarlo v0.2.2
[74087812] Random123 v1.3.1
[fb686558] RandomExtensions v0.4.3
[e6cf234a] RandomNumbers v1.4.0
[c84ed2f1] Ratios v0.4.0
[3cdcf5f2] RecipesBase v1.1.1
[01d81517] RecipesPipeline v0.3.2
[731186ca] RecursiveArrayTools v2.11.4
[f2c3362d] RecursiveFactorization v0.1.12
[189a3867] Reexport v1.0.0
[ae029012] Requires v1.1.3
[ae5879a3] ResettableStacks v1.1.0
[37e2e3b7] ReverseDiff v1.9.0
[79098fc4] Rmath v0.7.0
[47965b36] RootedTrees v1.0.0
[7e49a35a] RuntimeGeneratedFunctions v0.5.2
[476501e8] SLEEFPirates v0.6.20
[1bc83da4] SafeTestsets v0.0.1
[0bca4576] SciMLBase v1.13.4
[30cb0354] SciMLTutorials v0.9.0
[6c6a2e73] Scratch v1.0.3
[efcf1570] Setfield v0.7.0
[992d4aef] Showoff v1.0.3
[699a6c99] SimpleTraits v0.9.3
[ed01d8cd] Sobol v1.5.0
[b85f4697] SoftGlobalScope v1.1.0
[a2af1166] SortingAlgorithms v1.0.0
[47a9eef4] SparseDiffTools v1.13.2
[276daf66] SpecialFunctions v1.4.1
[860ef19b] StableRNGs v1.0.0
[aedffcd0] Static v0.2.4
[90137ffa] StaticArrays v1.2.0
[82ae8749] StatsAPI v1.0.0
[2913bbd2] StatsBase v0.33.8
[4c63d2b9] StatsFuns v0.9.8
[789caeaf] StochasticDiffEq v6.34.1
[7792a7ef] StrideArraysCore v0.1.11
[09ab397b] StructArrays v0.5.1
[d1185830] SymbolicUtils v0.11.2
[0c5d862f] Symbolics v0.1.25
[3783bdb8] TableTraits v1.0.1
[bd369af6] Tables v1.4.2
[5d786b92] TerminalLoggers v0.1.3
[a759f4b9] TimerOutputs v0.5.9
[0796e94c] Tokenize v0.5.16
[3bb67fe8] TranscodingStreams v0.9.5
[592b5752] Trapz v2.0.2
[a2a6695c] TreeViews v0.3.0
[5c2747f8] URIs v1.3.0
[3a884ed6] UnPack v1.0.2
[1986cc42] Unitful v1.7.0
[3d5dd08c] VectorizationBase v0.20.11
[81def892] VersionParsing v1.2.0
[19fa3120] VertexSafeGraphs v0.1.2
[44d3d7a6] Weave v0.10.8
[efce3f68] WoodburyMatrices v0.5.3
[ddb6d928] YAML v0.4.6
[c2297ded] ZMQ v1.2.1
[a5390f91] ZipFile v0.9.3
[e88e6eb3] Zygote v0.6.11
[700de1a5] ZygoteRules v0.2.1
[6e34b625] Bzip2_jll v1.0.6+5
[83423d85] Cairo_jll v1.16.0+6
[3bed1096] Cuba_jll v4.2.1+0
[7bc98958] Cubature_jll v1.0.4+0
[5ae413db] EarCut_jll v2.1.5+1
[2e619515] Expat_jll v2.2.10+0
[b22a6f82] FFMPEG_jll v4.3.1+4
[f5851436] FFTW_jll v3.3.9+7
[a3f928ae] Fontconfig_jll v2.13.1+14
[d7e528f0] FreeType2_jll v2.10.1+5
[559328eb] FriBidi_jll v1.0.5+6
[0656b61e] GLFW_jll v3.3.4+0
[d2c73de3] GR_jll v0.57.2+0
[78b55507] Gettext_jll v0.21.0+0
[7746bdde] Glib_jll v2.68.1+0
[e33a78d0] Hwloc_jll v2.4.1+0
[1d5cc7b8] IntelOpenMP_jll v2018.0.3+2
[aacddb02] JpegTurbo_jll v2.0.1+3
[c1c5ebd0] LAME_jll v3.100.0+3
[dd4b983a] LZO_jll v2.10.1+0
[dd192d2f] LibVPX_jll v1.9.0+1
[e9f186c6] Libffi_jll v3.2.2+0
[d4300ac3] Libgcrypt_jll v1.8.7+0
[7e76a0d4] Libglvnd_jll v1.3.0+3
[94ce4f54] Libiconv_jll v1.16.1+0
[4b2f31a3] Libmount_jll v2.35.0+0
[89763e89] Libtiff_jll v4.1.0+2
[38a345b3] Libuuid_jll v2.36.0+0
[856f044c] MKL_jll v2021.1.1+1
[e7412a2a] Ogg_jll v1.3.4+2
[458c3c95] OpenSSL_jll v1.1.1+6
[efe28fd5] OpenSpecFun_jll v0.5.4+0
[91d4177d] Opus_jll v1.3.1+3
[2f80f16e] PCRE_jll v8.44.0+0
[30392449] Pixman_jll v0.40.1+0
[ea2cea3b] Qt5Base_jll v5.15.2+0
[f50d1b31] Rmath_jll v0.3.0+0
[fb77eaff] Sundials_jll v5.2.0+1
[a2964d1f] Wayland_jll v1.17.0+4
[2381bf8a] Wayland_protocols_jll v1.18.0+4
[02c8fc9c] XML2_jll v2.9.12+0
[aed1982a] XSLT_jll v1.1.34+0
[4f6342f7] Xorg_libX11_jll v1.6.9+4
[0c0b7dd1] Xorg_libXau_jll v1.0.9+4
[935fb764] Xorg_libXcursor_jll v1.2.0+4
[a3789734] Xorg_libXdmcp_jll v1.1.3+4
[1082639a] Xorg_libXext_jll v1.3.4+4
[d091e8ba] Xorg_libXfixes_jll v5.0.3+4
[a51aa0fd] Xorg_libXi_jll v1.7.10+4
[d1454406] Xorg_libXinerama_jll v1.1.4+4
[ec84b674] Xorg_libXrandr_jll v1.5.2+4
[ea2f1a96] Xorg_libXrender_jll v0.9.10+4
[c7cfdc94] Xorg_libxcb_jll v1.13.0+3
[cc61e674] Xorg_libxkbfile_jll v1.1.0+4
[12413925] Xorg_xcb_util_image_jll v0.4.0+1
[2def613f] Xorg_xcb_util_jll v0.4.0+1
[975044d2] Xorg_xcb_util_keysyms_jll v0.4.0+1
[0d47668e] Xorg_xcb_util_renderutil_jll v0.3.9+1
[c22f9ab0] Xorg_xcb_util_wm_jll v0.4.1+1
[35661453] Xorg_xkbcomp_jll v1.4.2+4
[33bec58e] Xorg_xkeyboard_config_jll v2.27.0+4
[c5fb5394] Xorg_xtrans_jll v1.4.0+3
[8f1865be] ZeroMQ_jll v4.3.2+6
[3161d3a3] Zstd_jll v1.5.0+0
[0ac62f75] libass_jll v0.14.0+4
[f638f0a6] libfdk_aac_jll v0.1.6+4
[b53b4c65] libpng_jll v1.6.38+0
[a9144af2] libsodium_jll v1.0.20+0
[f27f6e37] libvorbis_jll v1.3.6+6
[1270edf5] x264_jll v2020.7.14+2
[dfaa095f] x265_jll v3.0.0+3
[d8fb68d0] xkbcommon_jll v0.9.1+5
[56f22d72] Artifacts
[2a0f44e3] Base64
[8bb1440f] DelimitedFiles
[8ba89e20] Distributed
[7b1f6079] FileWatching
[9fa8497b] Future
[b77e0a4c] InteractiveUtils
[4af54fe1] LazyArtifacts
[b27032c2] LibCURL
[76f85450] LibGit2
[8f399da3] Libdl
[37e2e46d] LinearAlgebra
[56ddb016] Logging
[d6f4376e] Markdown
[ca575930] NetworkOptions
[44cfe95a] Pkg
[de0858da] Printf
[9abbd945] Profile
[3fa0cd96] REPL
[9a3f8284] Random
[ea8e919c] SHA
[9e88b42a] Serialization
[1a1011a3] SharedArrays
[6462fe0b] Sockets
[2f01184e] SparseArrays
[10745b16] Statistics
[4607b0f0] SuiteSparse
[fa267f1f] TOML
[a4e569a6] Tar
[8dfed614] Test
[cf7118a7] UUIDs
[4ec0a83e] Unicode
[e66e0078] CompilerSupportLibraries_jll
[deac9b47] LibCURL_jll
[29816b5a] LibSSH2_jll
[c8ffd9c3] MbedTLS_jll
[14a3606d] MozillaCACerts_jll
[4536629a] OpenBLAS_jll
[bea87d4a] SuiteSparse_jll
[83775a58] Zlib_jll
[8e850ede] nghttp2_jll
[3f19e933] p7zip_jll

format: plain


text:
# Gravitational potential at the center of a uniform sphere

Feynman in second Messenger Lecture said that potential at the center of ball with small radius $a$ is equal to average potential on surface of ball minus $G$ times mass inside the ball divided by $2a$. You can see it on the picture. I don't understand that. I had an idea that it is from Taylor series of potential at the centre, but I can't figure it out. Can you help me?

• Welcome to physics.SE! Do you know the shell theorem, and do you know how to prove from it that the field inside a uniform sphere is $g=g_0r/a$, where $g_0$ is the field at the surface? – Ben Crowell Feb 19 '18 at 18:25
• I know shell theorem. Second I maybe get by substituting m=V.density and V=4/3 pi r^3. What do you mean? – Mat Fy Feb 19 '18 at 19:47
• Can you explain me second thing? Thank you for comment. – Mat Fy Feb 20 '18 at 6:11

Feynman's talking about a ball, which means that he is talking about a solid sphere, with uniform density, which I shall call $\rho$. You can apply Gauss's law for gravity to then calculate the potential.

Gauss's law states that:

$$\int FdA = -4\pi GM$$

where F is the g-field, A is a surface area and M is the mass enclosed by our Gaussian surface.

Let's say that our ball has radius $a$. We can imagine a Gaussian sphere, of radius $r < a$ which encloses some mass $\rho V$ where V is the volume of our Gaussian sphere (we can use multiplication instead of integration here because it's a uniform sphere). Call A the Gaussian area and by symmetry, $\int FdA = FA$, so that we can write, for our Gaussian sphere:

$$F = \frac {-4πG\rho V}{A}$$

For our Gaussian sphere, $V = \frac{4}{3}\pi r^3$ and $A = 4\pi r^2$ so that:

$$\frac VA = \frac r3$$

$$F = \frac{-4πG\rho r}{3}$$

This tells us the field at any point in the sphere. Replacing r with a tells us the field at the surface of the sphere.

Anyway, we want the potential at the surface of the sphere, which is the negative work done per unit mass when moving from the center of the sphere to its surface i.e from $r=0$ to $r=a$ (assuming that our reference point is set at the center of the sphere, so that the center has zero potential energy). This is the same thing as the integral of the field along this distance, since the field is the force per unit mass. Thus the potential at the surface of the sphere is:

$$\int_0^a Fdr =\int_0^a \frac{4πG\rho r dr}{3}$$

which is:

$$\frac{4πG\rho a^2}{6}$$

The volume of the entire sphere is $V = \frac 43 \pi a^3$, so that $\frac {4\pi a^2}{6} = \frac {V}{2a}$. Since $M$, the mass enclosed by the sphere is $\rho V$, where this time V is the volume of the entire sphere:

$$\frac{4πG\rho a^2}{6} = \frac{G\rho V}{2a} = \frac {GM}{2a}$$

So the work done in moving from the center to the surface of a ball is given by $\frac {GM}{2a}$, provided the ball is uniform. To get Feynman's equation just change the reference point for the potential energy, so that the potential energy at the center is no longer $0$.

Edit:

The above result holds only if the field due to matter outside the sphere is the same at the center and on the surface. If not the work done will depend on where on the surface you move to, since every point on the surface is not equidistant from external matter. If the ball is small enough, a criterion which Feynman mentions, then there's not much difference between different points on the surface or at the center, so you can approximate the field as constant on the ball.

• The work done in moving from the center to the surface of a ball is GM/2a only if potential energy at the center is zero. If you change field, then potential energy and so work will be different. So result is only approximate. – Mat Fy Jul 20 '18 at 14:24
• Ah yeah, Feynman does mention that this holds only if the ball is small enough, presumably so that the field due to external factors at the surface is the same as that at the center. I forgot to mention that. Thanks, I'll add it. – Albertrichard Jul 21 '18 at 14:45

Feynman's answer does not refer to a sphere filled with uniform mass density. It is a completely general statement about the potential at the center of a (small enough) sphere given an arbitrary mass distribution. It is a different way of stating Newton's law of gravitation. It can be derived from the differential formulation of Newton's law (ΔΦ = 4πG ρ) with the help of the full-space Green's function G(r,r')=-1/(4 π |r-r'|).

• Could you be more concrete, please? – Mat Fy Jan 30 '19 at 16:07


text:
Lemma 15.74.6. Let $R$ be a ring. Let $K^\bullet$ be a bounded complex of perfect $R$-modules. Then $K^\bullet$ is a perfect complex.

Proof. Follows by induction on the length of the finite complex: use Lemma 15.74.4 and the stupid truncations. $\square$

There are also:

• 6 comment(s) on Section 15.74: Perfect complexes

In your comment you can use Markdown and LaTeX style mathematics (enclose it like $\pi$). A preview option is available if you wish to see how it works out (just click on the eye in the toolbar).


text:
## $L^p$ restriction of eigenfunctions to random Cantor-type sets

Series
Analysis Seminar
Time
Wednesday, September 26, 2018 - 1:55pm for 1 hour (actually 50 minutes)
Location
Skiles 005
Speaker
Suresh Eswarathasan – Cardiff University
Organizer
Shahaf Nitzan
Abstract: Let $(M,g)$ be a compact Riemannian n-manifold without boundary. Consider the corresponding $L^2$-normalized Laplace-Beltrami eigenfunctions. Eigenfunctions of this type arise in physics as modes of periodic vibration of drums and membranes. They also represent stationary states of a free quantum particle on a Riemannian manifold. In the first part of the lecture, I will give a survey of results which demonstrate how the geometry of $M$ affects the behaviour of these special functions, particularly their “size” which can be quantified by estimating $L^p$ norms. In joint work with Malabika Pramanik (U. British Columbia), I will present in the second part of my lecture a result on the $L^p$ restriction of these eigenfunctions to random Cantor-type subsets of $M$. This, in some sense, is complementary to the smooth submanifold $L^p$ restriction results of Burq-Gérard-Tzetkov ’06 (and later work of other authors). Our method includes concentration inequalities from probability theory in addition to the analysis of singular Fourier integral operators on fractals.


text:
Best books on Set Theory?

1. Jan 10, 2008

andytoh

Any books that really stand out? Currently, I'm reading "Set Theory and Logic" by Stoll. I'm not interested in the axiomatic type of set theory, like Godel's theory and all those unreadable symboic proofs. I'm more interested in stuff like the axiom of choice proofs and such. Also, is there any research left to do in set theory or is it a fully exhausted field?

Last edited: Jan 10, 2008
2. Jan 11, 2008

Son Goku

I don't know about others, but I quite enjoyed Keith Devlin's "The Joy of Sets" (2nd Edition is much better). It gives a good justification of each of Zermelo's axioms and why they are there and has a very good explanation of the ordinals.
However the unique feature of this book is that it contains a good beginner's explanation of why $$2^{\aleph_{0}} = \aleph_{1}$$ is independant of ZFC and the attempts to resolve this by the addition of new axioms.

3. Jan 11, 2008

MathematicalPhysicist

well there are plenty of good intro books on set theory.
At my school we're using jech and hrbaceck intro to set theory, after reading this intro I think the next step is reading jech's set theory text, which is more advanced.

As for research in the field, if you search for handbook of set theory in google you'll find a page with the articles from the handbook, they address there their research in the field, so yes there's research in the field.


text:
## speed of light

$\Delta p \Delta x\geq \frac{h}{4\pi }$

Elyssa Smith 4B
Posts: 58
Joined: Wed Sep 18, 2019 12:22 am

### speed of light

Today in Lavelle's review, he said that we know the electron cannot be inside the nucleus of an atom because the calculated uncertainty in the velocity was greater than the speed of light. Why are we comparing it to the speed of light? Also, if it were smaller than 3x10^8, could it then be found inside the nucleus?

KeyaV1C
Posts: 78
Joined: Sat Aug 17, 2019 12:17 am
Been upvoted: 1 time

### Re: speed of light

Since no particles can travel faster than the speed of light, it makes no sense that a particle with mass (ie. electron) would be able to travel faster than the speed of light.

Aiden Metzner 2C
Posts: 72
Joined: Wed Sep 18, 2019 12:21 am

### Re: speed of light

Electrons are never found in the nucleus even if the uncertainty in velocity is less than the speed of light. That being said I don't think there is a case where this is ever possible because heisenbergs uncertainty principle would likely always show that an electrons velocity uncertainty is more then the speed of light if the uncertainty in position is confined in the nucleus.


text:
```Joseph Wright wrote:
> Will Robertson wrote:
>  > Another random thought: since users "shouldn't" be defining new macros
>> within the body of their documents, I think LaTeX3 should have something
>> like the equivalent of
>>   \AtBeginDocument{\catcode`\#=12}
>
> This of course would be a major break with LaTeX2e at the document
> level.  However, as everyone seems agreed that there is no real chance
> of compiling arbitrary LaTeX2e documents with a full LaTeX3 kernel, I'd
> say it is a design choice worth making.

But then I would say \catcode =12 to everything else, excluding only \
and maybe {} for grouping. Maybe even the meaning of % and \% could be
interchanged: % in normal Text and \% as comment?
Or is that too much changed?…

Cheers
Arno

```


text:
Here are the number of new confirmed COVID-19 cases per day, here in the province of Ontario:

In the great United States:

And in the entire world:

The plots tell the story. I don’t think that there is anything that I can add.

Dame Vera Lynn is no longer with us.

She was 103 years old.

In response to the COVID-19 pandemic, the Parliament of Hungary granted Viktor Orban’s government extraordinary powers to rule by decree. Opponents likened this to the infamous “enabling act” that marked the end of the Weimar Republic in 1933.

In response, Mr. Orban and his supporters argued that the government has no intention of abusing these powers, and that Mr. Orban should be judged by how he relinquishes these extraordinary powers once the crisis is over.

So here it is: Parliament just voted to end the state of emergency and requested the government to take the necessary steps to return to normalcy. So far so good. The corresponding act of Parliament is short and sweet, and to the point. The legal part of the text is only four terse clauses. In my translation:

1. The Legislature requests the government to end, in accordance with section 54, paragraph 3 of the Constitution, the state of emergency that was declared in government decree 40/2020 (March 11) (henceforth: the state of emergency).

2. The act of Parliament concerning measures against the coronavirus, Bill 12 (2020), is no longer in force.

3. (1) The present act – with the exceptions defined in paragraph (2) – shall be considered in force the day after it is proclaimed.

(2) Clauses 2 and 4 shall be considered in force when the state of emergency ends.

(3) The calendar day when clauses 2 and 4 shall be considered in force will be decided by the Prime Minister as soon as the date is known, by immediately proclaiming the date in the Hungarian Gazette.

4. Clause 2 of this bill, in accordance with section 29, paragraph 3; section 2, paragraph 1; section 24, paragraph 9, section 31, paragraphy 3; section 35, paragrapha 1; and section 54, paragraph 4, is considered fundamental.

All looks kosher, right?

But just in case you are still paying attention… there was also a second bill on the table. It has an unassuming title: “About the rules of transition concerning the end of the state of emergency and health preparedness”. Seems eminently reasonable, since under the state of emergency, the government has taken many steps, all of which must be reconciled with the system of laws and regulations.

And there is a lot. This supplementary bill is not four pages. It is 247 (two hundred forty seven) pages in PDF form. Yikes!

And on page 135, we read clause 97:

97. Changes to Bill 128 of 2011 concerning emergency preparedness and modifications of related laws

§339 Section 5 of Bill 128 of 2011 concerning emergency preparedness and modifications of related laws will be amended by the following subclause 24/A:

“24/A Emergency acts of government in case of a state of emergency concerning a human epidemic causing mass rates of infection”

§51/A. (1) The Government, in order to prevent, or mitigate the consequences of, a human epidemic that causes mass rates of infection, threatening lives and property, in order to protect the health and life of citizens of Hungary in a declared state of emergency – in addition to the emergency measures and rules described in subclauses 21-24 – in order to guarantee the lives, health, personal security, security of property, and legal protection of citizens and the stability of the national economy, may, by decree, suspend certain laws, deviate from existing laws and make other extraordinary decisions.

(2) The Government may practice the authority prescribed in paragraph (1) – to the necessary extent, in proportion with the desired outcome – in order to prevent, deal with, and eliminate an epidemic, and furthermore, to prevent or avoid such an outbreak.

There you have it. There are many other clauses concerning other laws, and most seem quite reasonable and appropriate considering the circumstances. But this one?

I highlighted the problem text in red. The law gives the government extraordinary powers. But, though the words seem reassuring on the surface, it offers no checks and balances concerning any justification of invoking these powers and the extent to which they are practiced. In the right hands, these are powers that can be used wisely in case of another pandemic or a resurgence of COVID-19, and I am sure Mr. Orban’s supporters believe that this is precisely how the government will use these powers (if at all). But opponents of Mr. Orban are concerned, not without grounds, about the lack of checks and balances and the very real possibility that these powers can go unchecked in the wrong hands.

In short, this is precisely what many feared: That though the state of emergency is lifted, a backdoor remains, a means for Mr. Orban’s government to have its cake and eat it, too, relinquish its extraordinary powers yet keeping them at the same time.

In the end, it of course all depends on how the law is put into practice. But for now, my conclusion is that the concerns of Mr. Orban’s opposition are not unfounded.

I just happened upon this 30-year old John Cleese sketch, where else? On John Cleese’s own Twitter feed, of course.

What can I say? Just perfect. As appropriate today as it was three decades ago, perhaps even more so.

In 1889, a story by Jules Verne (believed to have been written actually by his son, Michel Verne) was published in the American magazine Forum under the title, “In the Year 2889“.

In it, among other things, Verne envisions video conferencing.

Verne’s story was illustrated by George Roux, who is best known for his numerous illustrations for Verne’s science-fiction novels. I suspect that this particular picture was made in 1889 or 1890 (when Verne’s story, which appeared originally in English, was republished in France.)

I find this image mind-boggling. That 130 years ago, back in the 19th century, someone was able to envision… well, something that, for all intents and purposes, looks pretty much like what many of us are doing today.

Long overdue, but I just finished preparing the latest Maxima release, version 5.44.

I am always nervous when I do this. It is one thing to mess with my own projects, it is another thing to mess with a project that is the work of many people and contains code all the way back from the 1960s.

In case anyone doubted that modern birds are descendants of dinosaurs, here is a reminder: the shoebill.

These amazing creatures are apparently quite docile with humans, but eat baby crocodiles for lunch, which they kill by decapitating them.

They really look like survivors of the K-T asteroid impact. They are… I think they are beautiful.

Amidst all the tension that has been unleashed in the United States, there is this small ray of hope.

A black flight attendant on a Southwest flight initiated a conversation with a white passenger, who was reading the book White Fragility: Why It’s So Hard for White People to Talk About Racism, by Robin DiAngelo.

The white passenger’s remark, “It’s our fault. We have to start these conversations,” caught her by surprise. A short conversation followed. Then the big revelation: The unassuming gentleman happened to be Doug Parker, CEO of American Airlines.

I can already hear some of my friends objecting: “It’s not our fault!” Do not misconstrue Parker’s words (perhaps they weren’t even quoted verbatim.) He of course didn’t mean, I am sure, that every white person must bear personal responsibility for every vile act of racism that happens in America or elsewhere.

Rather, what I read into those words is an acknowledgement of a simple reality: In an unequal relationship, the dominant party has the power to make change for the better. In America, this means whites.

The fact that the CEO of a company as large as American Airlines recognizes this is, well, a ray of hope. As is the fact that he traveled, unassumingly, as an ordinary economy passenger on a competitor’s flight. As rising inequality between the super-wealthy and the stagnating middle class plagues Western societies, the US in particular, as disadvantaged minorities fall even further behind, it is nice to know that at least some folks in positions of power recognize that their wealth and status also come with a huge responsibility. Especially if the nice thoughts are also followed by deeds.

I don’t always like commercial publishers. Some of their textbooks are prohibitively expensive, yet often lacking in quality. (One persistent exception is Dover Publications, who published some of the best textbooks I own, as low-cost paperbacks.)

Last night, however, I was very pleasantly surprised by Springer, who made several hundred textbooks across a range of disciplines available for free, on account of COVID-19.

I did not get greedy. I didn’t download titles indiscriminately. But I did find several titles that are of interest to me, and I gladly took advantage of this opportunity.

Thank you, Springer.

Still playing with some COVID-19 maps, so here is another one: This one ranks countries by the number of COVID-19 cases per 1,000 square kilometers.

What’s the point, you might ask? Well, when there are lots of people confined to a small area, even a few cases can mean trouble; in contrast, when you have many cases but spread over half a continent, you may never come across a single infected person in your travels.

Again, the color legend remains a little whacky; it is logarithmic, but I don’t know how to convince this R package to display it nicely.

The numbers are pretty high. It goes without saying that densely populated microstates like Monaco win the contest, but then there is Qatar (3494), Belgium (1866), the Netherlands (1324), the UK (1051)… numbers that are way too high. For comparison, the USA is at 180, Hungary at 41, Russia at 20 (no surprise there), China at 9 (!), and so is Canada.

Maps that show the spread and mortality rate of COVID-19 are worrisome, to say the least.

So here is something slightly more encouraging: the day-to-day growth rate of COVID-19 infections, in the world, the United States, Canada, and the province of Ontario where I live.

What this plot shows is that the growth rate has been consistently and steadily decreasing (note that the vertical axis is logarithmic). Furthermore, the overall behavior of the worldwide, US, and Canadian plots is remarkable similar.

Of course this plot is also a reflection of the fact that the cumulative number of cases increases, so even if the number of new cases remains steady, the growth rate will indeed decrease. Nonetheless, it demonstrates that at least for now, the pandemic’s exponential spread has been arrested.  This can, of course, change for the worse quite dramatically if we lose our collective patience and start relaxing too soon.

I suppose the most important thing to know about COVID-19 is not so much the total number of cases or deaths, but the number of cases or deaths per million people. A number like 1,000 deaths is huge in a country of one million, a drop in the bucket in a country like India or China.

In that vein, I produced two maps using my rather rusty R programming skills. I am sure I could have done a better job labeling the legend, but there is only so much time I wanted to spend on this. The color bar legend is logarithmic between its two end values.

First, the cumulative number of cases per million people (May 22 data):

Next, deaths per million:

Not a pretty picture. It appears that the countries that are the most disproportionately affected (or the most likely to offer accurate reports?) are in Europe (including Russia), North and South America and the Middle East.

But now, here is perhaps the scariest map of all:

Yes, in many countries more than 10% of the confirmed cases result in death. The global average is 6.4 per 100 confirmed cases. The fact that COVID-19 is deadly in countries with underdeveloped health care systems is perhaps understandable, but the fact that the reddest parts of the map are places with the most sophisticated health care systems in the world is food for thought.

I am one of the maintainers of the Maxima computer algebra system. Maxima’s origins date back to the 1960s, when I was still in kindergarten. I feel very privileged that I can participate in the continuing development of one of the oldest continuously maintained software system in wide use.

It has been a while since I last dug deep into the core of the Maxima system. My LISP skills are admittedly a bit rusty. But a recent change to a core Maxima capability, its ability to create Taylor-series expansions of expressions, broke an important feature of Maxima’s tensor algebra packages, so it needed fixing.

The fix doesn’t amount to much, just a few lines of code:

It did take more than a few minutes though to find the right (I hope) way to implement this fix.

Even so, I had fun. This is the kind of programming that I really, really enjoy doing. Sadly, it’s not the kind of programming for which people usually pay you Big Bucks… Oh well. The fun alone was worth it.

One of the most fortunate moments in my life occurred in the fall of 2005, when I first bumped into John Moffat, a physicist from The Perimeter Institute in Waterloo, Ontario, Canada, when we both attended the first Pioneer Anomaly conference hosted by the International Space Science Institute in Bern, Switzerland.

This chance encounter turned into a 15-year collaboration and friendship. It was, to me, immensely beneficial: I learned a lot from John who, in his long professional career, has met nearly every one of the giants of 20th century physics, even as he made his own considerable contributions to diverse areas ranging from particle physics to gravitation.

In the past decade, John also wrote a few books for a general audience. His latest, The Shadow of the Black Hole, is about to be published; it can already be preordered on Amazon. In their reviews, Greg Landsberg (CERN), Michael Landry (LIGO Hanford) and Neil Cornish (eXtreme Gravity Institute) praise the book. As I was one of John’s early proofreaders, I figured I’ll add my own.

John began working on this manuscript shortly after the announcement by the LIGO project of the first unambiguous direct detection of gravitational waves from a distant cosmic event. This was a momentous discovery, opening a new chapter in the history of astronomy, while at the same time confirming a fundamental prediction of Einstein’s general relativity. Meanwhile, the physics world was waiting with bated breath for another result: the Event Horizon Telescope collaboration’s attempt to image, using a worldwide network of radio telescopes, either the supermassive black hole near the center of our own Milky Way, or the much larger supermassive black hole near the center of the nearby galaxy M87.

Bookended by these two historic discoveries, John’s narrative invites the reader on a journey to understand the nature of black holes, these most enigmatic objects in our universe. The adventure begins in 1784, when the Reverend John Michell, a Cambridge professor, speculated about stars so massive and compact that even light would not be able to escape from its surface. The story progresses to the 20th century, the prediction of black holes by general relativity, and the strange, often counterintuitive results that arise when our knowledge of thermodynamics and quantum physics is applied to these objects. After a brief detour into the realm of science-fiction, John’s account returns to the hard reality of observational science, as he explains how gravitational waves can be detected and how they fit into both the standard theory of gravitation and its proposed extensions or modifications. Finally, John moves on to discuss how the Event Horizon Telescope works and how it was able to create, for the very first time, an actual image of the black hole’s shadow, cast against the “light” (radio waves) from its accretion disk.

John’s writing is entertaining, informative, and a delight to follow as he accompanies the reader on this fantastic journey. True, I am not an unbiased critic. But don’t just take my word for it; read those reviews I mentioned at the beginning of this post, by preeminent physicists. In any case, I wholeheartedly recommend The Shadow of the Black Hole, along with John’s earlier books, to anyone with an interest in physics, especially the physics of black holes.

Heaven knows why I sometimes get confused by the simplest things.

In this case, the conversion between two commonly used cosmological coordinate systems: Comoving coordinates vs. coordinates that are, well, not comoving, in which cosmic expansion is ascribed to time dilation effects instead.

In the standard coordinates that are used to describe the homogeneous, isotropic universe of the Friedmann–Lemaître–Robertson–Walker (FLRW) metric, the metric is given by

$$ds^2=dt^2-a^2dR^2,$$

where $$a=a(t)$$ is a function of the time coordinate, and $$R$$ represents the triplet of spatial coordinates: e.g., $$dR^2=dx^2+dy^2+dz^2.$$

I want to transform this using $$R’=aR,$$ i.e., transform away the time-dependent coefficient in front of the spatial term in the metric. The confusion comes because for some reason, I always manage to convince myself that I also have to make the simultaneous replacement $$t’=a^{-1}dt.$$

I do not. This is nonsense. I just need to introduce $$dR’$$. The rest then presents itself automatically:

\begin{align*} R’&=aR,\\ dR&=d(a^{-1}R’)=-a^{-2}\dot{a}R’dt+a^{-1}dR’,\\ ds^2&=dt^2-a^2[-a^{-2}\dot{a}R’dt+a^{-1}dR’]^2\\ &=(1-a^{-2}\dot{a}^2{R’}^2)dt^2+2a^{-1}\dot{a}R’dtdR’-d{R’}^2\\ &=(1-H^2{R’}^2)dt^2+2HR’dtdR’-d{R’}^2, \end{align*}

where $$H=\dot{a}/a$$ as usual.

OK, now that I recorded this here in my blog for posterity, perhaps the next time I need it, I’ll remember where to find it. For instance, the next time I manage to stumble upon one of my old Quora answers that, for five and a half years, advertised my stupidity to the world by presenting an incorrect answer on this topic.

This, incidentally, would serve as a suitable coordinate system representing the reference frame of an observer at the origin. It also demonstrates that such an observer sees an apparent horizon, the cosmological horizon, given by $$1-H^2{R’}^2=0,$$, i.e., $$R’=H^{-1},$$ the distance characterized by the inverse of the Hubble parameter.

So here I am, reading about some trivial yet not-so-trivial probability distributions.

Let’s start with the uniform distribution. Easy-peasy, isn’t it: a random number, between 0 and 1, with an equal probability assigned to any value within this range.

So… what happens if I take two such random numbers and add them? Why, I get a random number between 0 and 2 of course. But the probability distribution will no longer be uniform. There are more ways to get a value in the vicinity of 1 than near 0 or 2.

And what happens if I add three such random numbers? Or four? And so on?

The statistics of this result are captured by the Irwin-Hall distribution, defined as

$$f_{\rm IH}(x,n)=\dfrac{1}{2(n-1)!}\sum\limits_{k=1}^n(-1)^k\begin{pmatrix}n\\k\end{pmatrix}(x-k)^{n-1}{\rm sgn}(x-k).$$

OK, so that’s what happens when we add these uniformly generated random values. What happens when we average them? This, in turn, is captured by the Bates distribution, which, unsurprisingly, is just the Irwin-Hall distribution, scaled by the factor $$n$$:

$$f_{\rm B}(x,n)=\dfrac{n}{2(n-1)!}\sum\limits_{k=1}^n(-1)^k\begin{pmatrix}n\\k\end{pmatrix}(nx-k)^{n-1}{\rm sgn}(nx-k).$$

For what it’s worth, here is the Maxima script to generate the Irwin-Hall plot:

fI(x,n):=1/2/(n-1)!*sum((-1)^k*n!/k!/(n-k)!*(x-k)^(n-1)*signum(x-k),k,0,n);
plot2d([fI(x,1),fI(x,2),fI(x,4),fI(x,8),fI(x,16)],[x,-2,18],[box,false],
[legend,"n=1","n=2","n=4","n=8","n=16"],[y,-0.1,1.1]);

And this one for the Bates plot:

fB(x,n):=n/2/(n-1)!*sum((-1)^k*n!/k!/(n-k)!*(n*x-k)^(n-1)*signum(n*x-k),k,0,n);
plot2d([fB(x,1),fB(x,2),fB(x,4),fB(x,8),fB(x,16)],[x,-0.1,1.1],[box,false],
[legend,"n=1","n=2","n=4","n=8","n=16"],[y,-0.1,5.9]);

Yes, I am still a little bit of a math geek at heart.

My lovely wife, Ildiko, woke up from a dream and asked: If you have a flower with 7 petals and two colors, how many ways can you color the petals of that flower?

Intriguing, isn’t it.

Such a flower shape obviously has rotational symmetry. Just because the flower is rotated by several times a seventh of a revolution, the resulting pattern should not be counted as distinct. So it is not simply calculating what number theorists call the $$n$$-tuple. It is something more subtle.

We can, of course, start counting the possibilities the brute force way. It’s not that difficult for a smaller number of petals, but it does get a little confusing at 6. At 7 petals, it is still something that can be done, but the use of paper-and-pencil is strongly recommended.

So what about the more general case? What if I have $$n$$ petals and $$k$$ colors?

Neither of us could easily deduce an answer, so I went to search the available online literature. For a while, other than finding some interesting posts about cyclic, or circular permutations, I was mostly unsuccessful. In fact, I began to wonder if this one was perhaps one of those embarrassing little problems in combinatorial mathematics that has no known solution and about which the literature remains strangely quiet.

But then I had another idea: By this time, we both calculated the sequence, 2, 3, 4, 6, 8, 14, 20, which is the number of ways flowers with 1, 2, …, 7 petals can be colored using only two colors. Surely, this sequence is known to Google?

Indeed it is. It turns out to be a well-known sequence in the online encyclopedia of integer sequences, A000031. Now I was getting somewhere! What was especially helpful is that the encyclopedia mentioned necklaces. So that’s what this problem set is called! Finding the Mathworld page on necklaces was now easy, along with the corresponding Wikipedia page. I also found an attempt, valiant though only half-successful if anyone is interested in my opinion, to explain the intuition behind this known result:

$$N_k(n)=\frac{1}{n}\sum_{d|n}\phi(d)k^{n/d},$$

where the summation is over all the divisors of $$n$$, and $$\phi(d)$$ is Euler’s totient function, the number of integers between $$1$$ and $$d$$ that are relative prime to $$d$$.

Evil stuff if you asked me. Much as I always liked mathematics, number theory was not my favorite.

In the case of odd primes, such as the number 7 that occurred in Ildiko’s dream, and only two colors, there is, however, a simplified form:

$$N_2(n)=\frac{2^{n-1}-1}{n}+2^{(n-1)/2}+1.$$

Substituting $$n=7$$, we indeed get 20.

Finally, a closely related sequence, A000029, characterizes necklaces that can be turned over, that is to say, the case where we do not count mirror images separately.

Oh, this was fun. It’s not like I didn’t have anything useful to do with my time, but it was nonetheless a delightful distraction. And a good thing to chat about while we were eating a wonderful lunch that Ildiko prepared today.

The other day, I ran across a question on Quora asking why Einstein didn’t support his country, Germany, during the Second World War. Thinking about this question reminded me of an old Star Trek episode and one of the root concepts (or, at least, my reading of it) of the Abrahamic family of religions.

In answering the question, I pointed out the difference between supporting a country vs. supporting a regime. I argued that Einstein, though not even a citizen of Germany at the time (he gave up German citizenship after Hitler’s rise to power in 1933, and became a naturalized US citizen in 1940), did, in fact, support his country of birth, precisely by the act of following his conscience and opposing the despotic, murderous Nazi regime.

And that takes me to the Star Trek episode Bread and Circuses from 1968. In this episode, the USS Enterprise encounters a planet governed by a regime not unlike the Roman Empire, but with 20th century technology, broadcasting gladiatorial matches by analog television. In due course, the crew of the Enterprise gets into trouble and link up with a group of rebellious Sun-worshippers. When at the end of the episode, after the conflict is resolved and the good guys prevail as usual, Spock expresses surprise over the fact that such a primitive religion could have survived on this planet into its modern era, Uhura corrects him by clarifying that they were, in fact, worshippers of the son of God. In other words, this planet’s version of early Christianity arrived two thousand years later than on the Earth.

Christianity borrows its creation mythology from Judaism, including the notion of the Garden of Eden and the Tree of Knowledge, the fruit of which let Adam and Eve understand the difference between good and evil. In my reading, this is what it really means when the Bible proclaims that humans are created in God’s image: that just like God, humans are free agents with a conscience, capable of acting independently, not robots blindly executing a predetermined divine script. They even have the capacity to act against God’s will.

Think about this, just what a revolutionary, what a deeply subversive concept this really is even today, never mind ancient times. The Book of Genesis is probably about 3,000 years old if not older. Egypt, in its third intermediate period, was ruled by pharaohs, seen as intermediaries between gods and ordinary people, whose words must be obeyed. Whether or not the Egyptian captivity happened (there do appear to be reasons to doubt), it’s no wonder Egypt’s rulers didn’t look kindly upon these pesky Jews and their subversive religion that claimed that it is more important to listen to your conscience than to blindly follow the orders of your divine ruler.

Despots can claim whatever they want: They can claim to represent the state, they can even claim to be the earthly representative of a divine power, like the pharaohs of old, but you have something over which they have no power: your conscience, which allows you to defy the will of any ruler, even God’s will, just as Adam and Eve have done back in the Garden of Eden.

And this is precisely what Einstein did when he lent his support, for instance, to Leo Szilard’s letter to Roosevelt that arguably launched the Manhattan project: Instead of slavishly following a despot claiming to represent the country of his birth, he listened to his conscience.

Just as things are beginning to look ever so slightly hopeful with infection rates at least stabilizing, conspiracy theorists are now having a field day.

In case you are wondering, we “know” that SARS-CoV-2 was manufactured (or at least released) by that Wuhan “bioweapons” lab, as its intended purpose was to weaken China’s strategic opponents, in particular the US military, so that they can mess with Taiwan. And even more blatantly, they did so using in part funding from the United States, according to information that was “just revealed”. And all in the service of some demonic Chinese plot to achieve some nefarious goal, such as the subjugation of the renegade province of Taiwan.

Like all good conspiracy theories, this nonsense is also based on a carefully picked set of selected facts. So much so that it reminds me of those old Soviet-era Radio Yerevan jokes.

Meanwhile in the real world…

1. We know that SARS-CoV-2 is “not a laboratory construct or a purposefully manipulated virus”. [Andersen, K.G., Rambaut, A., Lipkin, W.I. et al. The proximal origin of SARS-CoV-2. Nat Med 26, 450–452 (2020).]
2. The funding provided in part by the NIH was never a secret in the first place. The funding source is explicitly listed in research published by scientists from the Wuhan Institute of Virology. And while the research was indeed about SARS-related coronaviruses, its aim was (and remains) to help prevent or mitigate the illness, not cause it. In fact, some of the recipients of this grant went so far as to warn us that their “work provides new insights into the origin and evolution of SARS-CoV and highlights the necessity of preparedness for future emergence of SARS-like diseases.” [Hu B, Zeng L-P, Yang X-L, Ge X-Y, Zhang W, Li B, et al. (2017) Discovery of a rich gene pool of bat SARS-related coronaviruses provides new insights into the origin of SARS coronavirus. PLoS Pathog 13(11): e1006698.]
3. Yes, the two nuclear powered US aircraft carriers usually patrolling the Pacific are presently in port because of COVID-19 cases on board. Yes, China is parading their toy carrier in the region. But an amphibious assault ship, the USS America, is also in the region, and in case “amphibious assault ship” sounds less impressive than “aircraft carrier”, it may help to note that the USS America is similar in size to China’s Lianoning and unlike its Chinese counterpart, it is a modern warship commissioned in 2012, not something built upon the hull of an unfinished 1980s Soviet carrier that was once sold for scrap.

None of this matters to conspiracy theorists, of course. If they want to believe that it was all some evil Chinese plot, they will do so, damn the evidence. (Sadly, I think it is a safe bet that they get at least some help from Russian troll farms. Unfortunately Putin’s regime has not stopped spreading disinformation about health care, epidemics, pandemics, vaccines and now, the origins of COVID-19.)

So to all the conspiracy theory fans out there, please keep this in mind: Conspiracy theories serve a single purpose. They turn us: individuals or entire nations, against each other. They are the means to sow discord. Just when the world has to act in unison more than ever, they pull us apart. Conspiracy theories are weapons: dangerous in the wrong hands even when used unintentionally, but deadly if used purposefully.

Please… listen to the actual science. Don’t believe every piece of garbage you hear.

I haven’t blogged in two weeks. In my excuse, I was rather busy. The good kind of busy, that is, busy with paying work, busy with scientific research, not busy with illness or anything on that front.

Anyhow, even though I haven’t blogged, I’ve been keeping track of the numbers. And for the past few days, a ray of hope began to emerge.

To make a long story short, there are significant signs that mitigation measures are working. Here, this chart shows the doubling rate of COVID-19 infections worldwide:

Infection rates doubled every 5-6 days back in late March; now, the doubling rate is over 15 days and rapidly rising.

Perhaps the world data are manipulated. But then, here is the doubling rate for Canada. The data are much more noisy (the population is much smaller, so this is to be expected) but the similar trend is unmistakable:

But there is another sign that things just might be working. I’ve been following US data in more detail, and lately, the simple SIR model’s predictions began to match the data rather well. If the model is to be believed, we may not be out of the woods quite yet, but we may be surprisingly close:

What a difference a few weeks make: Back in late March, the same model predicted catastrophic numbers. Now, it seems to tell us that we are mere weeks away from life gradually beginning to return to normal.

I dare not believe it just yet, but it is a ray of hope.


text:
Consider the following limit We can simplify this limit by rewriting it as an expression without absolute values as follows
$\lim_{x\to 3}$
We can then cancel off a common factor in the numerator and denominator, thus simplifying our limit to
$\lim_{x\to 3}$
We can then evaluate the limit directly and find that its value is