drosophila embryo protein gradients

main.lua

--[[
	cartoon model of drosophila embryo protein gradients

	1,2 are basic anterior-posterior gradient
	1) red protein ~ hunchback
	2) green protein ~ caudal

	3,4,5 are 'gap genes'
	3) blue protein ~ kruppel

	6 is a pair rule gene, like even-skipped
	it only produces 3 stripes, but it works as a proof of concept


	run this with love2d   https://love2d.org/
]]

io.stdout:setvbuf("no")

width = 720
height = 720

love.window.setMode( width, height, {vsync = true} )

-- nr of cells
len = 100
-- nr of proteins
num = 6


colors = {
	{0.8,0.1,0.1},
	{0.1,0.8,0.1},
	{0.1,0.1,0.9},
	{0.8,0.8,0.1},
	{0.8,0.1,0.8},
	{0.1,0.8,0.8},
	{0.8,0.5,0.1},
}


function love.load()
	math.randomseed(os.time())
	local font = love.graphics.newFont( 17 )
	love.graphics.setFont(font)

	love.graphics.setLineWidth(1.0)
	love.graphics.setLineStyle("smooth")

	-- set up initial concentrations etc
	proteins = {}
	new = {}
	for i=1, num do
		proteins[i] = {}
		new[i] = {}
		for j=1,len do
			proteins[i][j] = 0
			new[i][j] = 0
		end
	end

	interactions = {}
	for i=1, num do
		interactions[i] = {}
		for j=1, num do
			interactions[i][j] = 0
		end
	end

	-- interaction terms
	-- positive is activator
	-- negative is repressor

	interactions[1][3] = 100
	interactions[2][3] = 100
	interactions[4][3] = -2
	interactions[5][3] = -2


	interactions[1][4] = 9
	interactions[2][4] = -28
	interactions[3][4] = 10

	interactions[2][5] = 12
	interactions[1][5] = -32
	interactions[3][5] = 11

	interactions[1][6] = -3
	interactions[2][6] = -3
	interactions[3][6] = 500
	interactions[4][6] = -3
	interactions[5][6] = -3


	-- diffusion and decay constants
	diffusion = {}
	decay = {}
	for i=1, num do
		diffusion[i] = 0.1
		decay[i] = 0.0005  
	end

	decay[3] = 0.005
	decay[4] = 0.003
	decay[5] = 0.004
	decay[6] = 0.008

	diffusion[4] = 0.03
	diffusion[5] = 0.03
	diffusion[6] = 0.05


	-- iterations per frame
	maxiter = 10
end

function love.update()
	mouseX,mouseY = love.mouse.getPosition()

	dt = 1

	for iter = 1,maxiter do
		for i,v in ipairs(proteins) do
			for j,b in ipairs(v) do
				local l = v[j-1]
				local r = v[j+1]
				local c = v[j]
				local n = 0

				-- do diffusion
				if l and r then
					n = n + diffusion[i]*(l+r-2*c)
				elseif l then
					n = n + diffusion[i]*(l-c)
				elseif r then
					n = n + diffusion[i]*(r-c)
				end

				--[[
					transcription activators / repressors
					note that the basal expression level is 1!
					the interactions are multiplicative,
					so if you have two activators, they need to be
					both present!
				]]
				local activity = 1
				for k = 1,num do
					local interaction = interactions[k][i]

					local a = interaction*proteins[k][j]

					-- these are Hill equation responses, with exponent 2  (kind of arbitrary)
					-- https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)#Regulation_of_gene_transcription
					if i~= k then
						if interaction > 0.0 then
							activity = activity*(a*a / (1+a*a))
						else
							activity = activity*(1 / (1+a*a))
						end
					end
				end

				if i > 2 then
					-- add protein based on transcription
					n = n + 0.01*activity
					-- this cubic term is simulating autokatalytic activity
					-- it makes the protein levels bistable
					-- in this case it also conventiently squashes the levels to [0-1]
					n = n + 0.05*(1-c)*c*(c-0.5)
					-- breakdown
					n = n - c * decay[i]
				else
					-- simpler model for the first two, basic gradients 
					n = n - c * decay[i]
				end

				-- add some noise
				--n = n + 0.0003*(math.random()-0.5)

				-- fwd euler update
				new[i][j] = c + dt*n

				if new[i][j] < 0 then
					new[i][j] = 0
				end
			end
		end
		-- polarity cheat
		-- the first and last cell produce 1st and 2nd protein resp
		new[1][1] = 1
		new[2][len] = 1

		-- swap buffer
		proteins, new = new, proteins
	end
end


function love.draw()
	love.graphics.setBackgroundColor(0.14,0.14,0.14)

	love.graphics.setColor(1, 1, 1)

	local sx = 6
	local sy = 200
	local x = 16
	local y = 300
	for i,v in ipairs(proteins) do
		if colors[i] then
			love.graphics.setColor(colors[i])
		end
		for j,b in ipairs(v) do
			if j < len then
				love.graphics.line(x + sx*j, y - sy*v[j], x + sx*(j+1), y - sy*v[j+1])
			end
		end
	end
end

function love.keypressed(key)
	if key == "r" then
		--reset
		for i=1, num do
			for j=1,len do
				proteins[i][j] = 0
			end
		end
	end

	if key == 'escape' then
		love.event.quit()
	end
end