Зуб комодского дракона

2024-11-03 в 9:00, admin, рубрики: python, визуализация, драконы, зубы, моделирование

Введение

Оценивать в миллиграммах содержание в ткани того или иного микроэлемента будем по формуле , где — плотность, V — объем.

Чтобы рассчитать объем, нужно знать форму объекта. А лучше — уравнение поверхности, его ограничивающей. С этого и начнем. Вообще задача моделирования поверхности интересна сама по себе, поэтому решим ее в нескольких вариантах, переходя от более простого к более точному.

Геометрия зуба

В качестве экспериментальных данных возьмем фотографии и описание из недавнего исследования, опубликованного в журнале Nature Ecology & Evolution [3]:

Зубы варана V. komodoensis. Источник изображения LeBlanc et al. Nat Ecol Evol 8, 1711–1722 (2024) по лицензии СС BY 4.0 — Зубы варана V. komodoensis. Источник изображения LeBlanc *et al.* *Nat Ecol Evol* 8, 1711–1722 (2024) по лицензии СС BY 4.0

На изображении видно, что зубы варана имеют загнутую и слегка сплющенную с боков форму. Длина зубов 2.5–3 см (что также указано в [1]). Зафиксируем для дальнейших расчетов высоту зуба равной 25 мм, глубину 10 мм, ширину 4 мм.

Зуб варана v.1 (конус)

Объем конического зуба

В качестве первого приближения возьмем конус с эллиптическим основанием.

Объем конуса, как известно со школы, $V=frac{1}{3} h S$

где S — площадь основания, h — высота.
В нашем случае h = 25 мм

В случае эллипса ,

где a и b — большая и малая полуоси эллипса, лежащего в основании.
В нашем случае a = 5 мм, b = 2 мм.

Таким образом, $V=frac{1}{3} h S=frac{1}{3} pi h a b$

Подставим численные значения: $V=frac{1}{3} pi cdot 25 cdot 2 cdot 5=261.7$

Итак объем зуба V = 261.7 мм³. Это наша первая грубая оценка.

Уравнение поверхности и визуализация

Способ 1 (быстрый)

Изобразим конус в логике аналитической геометрии, т.е. опираясь на уравнение его поверхности.

Уравнение конуса запишется как $frac{x^2}{a^2} + frac{y^2}{b^2}-frac{z^2}{h^2}=0$

Или, после подстановки значений a, b и h: $frac{x^2}{4} + frac{y^2}{25}-frac{z^2}{652}=0$

Чтобы задать конечную высоту зубу (сам по себе конус фигура бесконечная), наложим ограничение на диапазон изменения z:

Для быстрой визуализации можем воспользоваться, например, сервисом desmos.com. Вводим полученное уравнение как есть и наблюдаем изображение конуса.

Способ 2 (с полным контролем)

Изобразим конус при помощи Python и Matplotlib.

В данном случае будет удобнее взять уравнение конуса в параметрическом виде:

$begin{cases} z=t \ x=frac{a}{h} cdot t cdot cos(phi) \ y=frac{a}{h} cdot t cdot sin(phi) end{cases}$

где параметры изменяются в следующих диапазонах: ,

Код:

import numpy as np 
import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import Axes3D 

# Параметры конуса 
a, b, h = 2, 5, 25

# Координаты точек на поверхности конуса 
phi = np.linspace(0, 2 * np.pi, 30) 
tz = np.linspace(-h, 0, 20) 
phi, tz = np.meshgrid(phi, tz)

# Параметрическое уравнение конуса 
z = tz 
x = (a / h) * tz * np.cos(phi) 
y = (b / h) * tz * np.sin(phi)

# Создание графика 
ax = plt.figure().add_subplot(111, projection='3d') 
ax.view_init(elev=20, azim=-40, roll=0)
ax.plot_surface(x, y, z, color='w', edgecolor='royalblue',             	rstride=2, cstride=2, alpha=0.5) ax.set_aspect('equal')

# Диапазоны и настройки графика 
x_lims = [-2a, 2a] 
y_lims = [-1.5b, 1.5b] 
z_lims = [-h, 0]

ax.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
       xlabel='X', ylabel='Y', zlabel='Z',
       xticks=np.round(np.linspace(*x_lims, 5)))

plt.show()

В результате получим следующий график:

Однако поверхность зуба все же изогнутая, а не линейная как у конуса, поэтому стоит добавить кривизны. Понизим степень в третьем слагаемом и получим уравнение параболоида.

Зуб варана v.2 (параболоид)

На сей раз начнем с уравнения поверхности.

Уравнение поверхности и визуализация

Уравнение параболоида имеет вид $frac{x^2}{a^2} + frac{y^2}{b^2}+frac{z}{h}=0$

При подстановке заданных параметров получаем следующее уравнение боковой поверхности зуба:

$frac{x^2}{4} + frac{y^2}{25}+frac{z}{25}=0$

Как и прежде, ограничим диапазон изменения z:

Вводим полученное уравнение в сервисе desmos.com или аналогичном и наблюдаем соответствующий график параболоида.

Объем параболического зуба

Найдем объем зуба параболической формы как тройной интеграл по области пространства, ограниченной поверхностью параболоида и плоскостью .

В общем случае в выражении через повторные интегралы это выглядит так

$V=intlimits_{x_1}^{x_2}dx intlimits_{y_1(x)}^{y_2(x)}dy intlimits_{z_1(x,y)}^{z_2(x,y)}dz$

Чтобы удобно расставить пределы интегрирования выполним с помощью Matplotlib дополнительный чертеж.

Код для дополнительного чертежа

import matplotlib.pyplot as plt 
from mpl_toolkits.mplot3d import axes3d 
import numpy as np 

a, b, h = 2, 5, 25

# Создание поля для графика 
fig = plt.figure(constrained_layout=True, figsize=(7,5)) 
spec = fig.add_gridspec(ncols=2, nrows=1) 
spec_right = spec[1].subgridspec(ncols=1, nrows=5)

ax1 = fig.add_subplot(spec[0, 0], projection='3d') 
ax2 = fig.add_subplot(spec_right[1:4])

# ---------- График-1: Параболоид в 3D -----------

ax1.view_init(elev=20, azim=-40, roll=0)

# диапазоны по осям и сетка 
x_lims = [-1.5a, 1.5a] 
y_lims = [-1.5b, 1.5b] 
z_lims = [-1.0*h, 0]

X, Y = np.linspace(*x_lims, 50), np.linspace(*y_lims, 50) 
X, Y = np.meshgrid(X, Y)

def z_func(X, Y):
  return np.where((X2 / a2 + Y2 / b2) <= 1,
                  -h * (X2 / a2 + Y2 / b2) ,
                  -h)

Z = z_func(X,Y)

# Рисуем 3D поверхность 
ax1.plot_surface(X, Y, Z, edgecolor='royalblue', lw=0.5, 
                 rstride=3, cstride=3,
                 alpha=0.3) 
ax1.set_aspect('equal')

# Проекции на координатные плоскости 
ax1.contour(X, Y, Z, levels=[-h], zdir='z', offset=z_lims[0],
            colors=['black'])

# Подписи к графику 
ax1.text(x_lims[0]-1, y_lims[0], -1, r'', (0, 1, 0)) 
ax1.text(x_lims[0]-1, y_lims[0], -h, r'', (0, 1, 0))

ax1.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
        xlabel='X', ylabel='Y', zlabel='Z',
        xticks=np.linspace(*x_lims, 4))

# ------ График-2: Проекция на плоскость основания ------

# Эллипс 
phi = np.linspace(0, 2*np.pi, 50)

X = a * np.cos(phi) Y = b * np.sin(phi)

ax2.plot(X, Y, '--')

# Четвертинка 
X2 = np.linspace(0, a, 50) 
Y2 = b * np.sqrt(1 — X22 / a2)

ax2.plot(X2, Y2, color='orange') 
ax2.plot(X2, Y2*0, color='orange') 
ax2.fill_between(X2, Y2, 0, color='blue', alpha=.1)

# Аннотации 
ax2.text(0,-0.6, r'', rotation='vertical') 
ax2.text(a-0.2,-0.6, r'', rotation='vertical')
ax2.text(a + 0.2, 0, r'') ax2.text(a + 0.2, b, r'')

# Настройки отображения 
ax2.grid() 
ax2.set(xlabel=r'X', ylabel=r'Y')

plt.show()

Поскольку фигура симметричная, достаточно получить объем одной четверти, а потом умножить на 4. Это позволит выставить часть приделов равными 0, упростив тем самым выражение.

Таким образом для вычисления одной четверти объема:

интеграл по dz берем в пределах от поверхности (плоскости) основания
до поверхности параболоида $z_2(x,y)=-h (frac{x^2}{a^2}+frac{y^2}{b^2})$ ,
в интеграле по dy смотрим на проекцию параболоида на плоскость основания (получается эллипс) и берем в пределах
от прямой
до эллиптической границы основания $y_2(x)=bsqrt{1-frac{x^2}{a^2}}$ ,
в интеграле по dx смотрим на проекцию основания на ось координат Ox берем в пределах от точки до самой правой точки проекции

Объем всей фигуры тогда запишется как

$V=4intlimits_{0}^a dx intlimits_{0}^{bsqrt{1-frac{x^2}{a^2}}}dy intlimits_{-h}^{-h (frac{x^2}{a^2}+frac{y^2}{b^2})}dz$

Вычислим его с помощью функции tplquad() библиотеки SciPy:

from scipy import integrate 

a, b, h = 2, 5, 25

# Подынтегральная фунция и пределы интегрирования 
f = lambda x,y,z: 1 
x_lims = [0, a] 
y_lims = [0, lambda x: b * (1 — x*2 / a**2)0.5] 
z_lims = [-h, lambda x, y: -h * (x2 / a2 + y2 / b*2)]

# Тройной интеграл 
V, err = integrate.tplquad(f, *x_lims, *ylims, *zlims)

print(f'Объем = {V * 4:.2f}')

Объем зуба с параболическими обводами равен V=392.7 мм³ .

Продолжим совершенствование формы зуба и сделаем его более острым и изогнутым

Зуб варана v.3 (логарифм)

Уравнение поверхности и визуализация

Нам нужно «загнуть» центральную ось нашего параболического зуба. То есть, если посмотреть на рисунок ниже, — перейти от ситуации с левой верхней схемы к левой нижней.

Код рисунка

import sympy as sp
import numpy as np
from sympy.abc import x, y, z
from spb import *
from matplotlib.gridspec import GridSpec
from sympy.printing.latex import LatexPrinter


a, b, h = 2, 5, 25
la_print = LatexPrinter()  # рендеровщик LaTex формул

# Уравнение параболоида с логарифмическим сдвигом по "y"
left_part_eq = x**2/a**2 + (y + 2*sp.log(1-z))**2/b**2 + z/h
equations = [left_part_eq]

# Выразим "y" из уравнения
sol1, sol2 = sp.solve(equations, y, dict=True)

# Достанем программно логарифм из уравнения "загнутого" зуба
log_line = left_part_eq.args[2].args[1].args[0].args[1]
log_line

# Кривая смещения
eq_log = sp.Eq(y, -log_line)

# Уравнение параболической поверхности зуба
parab_left = x**2 / a**2 + y**2 / b**2 + z / h 
z_sol_parab = sp.solve([parab_left], z)

# Проекция параболоида
eq_parab = sp.Eq(z, z_sol_parab[z].subs({x:0}))

# Параметры отрисовки
params = {'aspect':(1,1), 'legend':False, 'show':False}
params_proj = {'aspect':(1,1), 'legend':False, 'show':False}
ranges_par = [(y, -4*b, 2*b), (z, -h, 0)]
ranges_log = [(y, -4*b, 2*b), (z, -h, 0)]
x_range = (x, -4*a, 4*a)
log_style = {"linestyles": "dashed", "linewidths": 1}
bacolor = 'dodgerblue'

# -------- График-1: Параболический зуб в профиль --------

# -- Аннотации
shift_arrows = []
for z_val in np.linspace(-0.2*h, -0.8*h, 3):
  y_val = -log_line.subs({z:z_val}).n()
  arrow = {'text': '', 'xy': (y_val, z_val), 
           'xytext': (0, z_val),
           'arrowprops': dict(arrowstyle="->")}
  shift_arrows.append(arrow)

# -- Парабола, вертикальное сечение зуба
p1 = plot_implicit(z < z_sol_parab[z].subs({x:0}),
                   *ranges_par, **params, 
                   xlabel='y', ylabel='z')

# -- Кривая куда смещаем, искривляя, ось
p1_logos = plot_implicit(eq_log, *ranges_par, color='red',
                         **params, annotations=shift_arrows, 
                         rendering_kw=log_style,)

# -- Кастомная легенда
legend_xy = (ranges_par[0][1], -2)
pretty_log_leg = la_print.doprint(sp.Eq(y,-log_line)).replace('log', 'ln')

text_legend = {'text': f'n${pretty_log_leg}$',
               'xy': (0, 0), 'xytext': legend_xy, 
               'fontsize': 8,
               'bbox': dict(boxstyle="square", fc="w")}

# -- Текущая прямая ось зуба-парабалода
p1_os = plot_implicit(sp.Eq(y, 0), *ranges_par, **params, 
                      color="white",
                      rendering_kw={"linestyles": "-.", 
                                    "linewidths": 4})

# -- График-2: Проекции на основание параболического зуба --

# -- Эллипс, проекция на основание
p1_base = plot_implicit(parab_left.subs({z:0}) <= 1,
                        (y, ranges_par[0][1], ranges_par[0][2]), x_range,
                        markers=[{'args': [0, 0, 'x'], 
                                  'color':'white'}],
                        #annotations=[text_para],
                        **params_proj, color=bacolor)

# -- Стрелка направление смещения основания в
new_cent = -log_line.subs({z:-h}).n()  # центр нового основания
shift_arrow_base = {'text': '', 'xy': (new_cent, 0), 'xytext': (0, 0),
                    'arrowprops': dict(arrowstyle="->")}
size_arrow_base = {'text': '', 'xy': (new_cent, 3), 'xytext': (0, 3),
                   'fontsize': 7,
                    'arrowprops': dict(arrowstyle="|-|")}
pretty_log = la_print.doprint(log_line).replace('log', 'ln')
text_log = {'text': f'${pretty_log}$',
             'xy': (0, 0), 'xytext': (new_cent, 5)}

# -- Эллипс, смещенная в новое положение проекция
p1_base_sh = plot_implicit(sp.Eq(left_part_eq.subs({z:-h}),0),
                           (y, ranges_par[0][1], ranges_par[0][2]), x_range,
                           annotations=[shift_arrow_base, size_arrow_base, text_log],
                           rendering_kw={"linestyles": "--", "linewidths": 1},
                           markers=[{'args': [new_cent, 0, 'x'], 'color':'red'}],
                           color='gray', **params_proj,
                           xlabel='y', ylabel='x')

# ------- График-3: Логарифмический зуб в профиль --------

# -- Вертикальное сечение изогнутого зуба
p2 = plot_implicit(sp.And(y >= sol1[y].subs({x:0}), y <= sol2[y].subs({x:0})),
                       *ranges_log, **params)

# -- Центральная ось изогнутого зуба
p2_os = plot_implicit(eq_log, *ranges_par, **params, color="white",
                      rendering_kw={"linestyles": "-.", "linewidths": 4})

# -- Кривая смещения (куда таки сместили, изогнув, параболоид)
p2_logos = plot_implicit(eq_log, *ranges_log, color='red',
                         **params, rendering_kw=log_style)

# -- График-4: Проекции на основание логарифмического зуба --

# -- Эллипс, проекция на основание
p2_base = plot_implicit(left_part_eq.subs({z:-h}) <= 0,
                        (y, ranges_log[0][1], ranges_log[0][2]), x_range,
                        markers=[{'args': [0, 0, 'x'], 'color':'gray'}],
                        #annotations=[text_loga],
                        **params_proj, color=bacolor)

# -- Эллипс, старая проекция, откуда смещались
p2_base_sh = plot_implicit(sp.Eq(parab_left.subs({z:-h}),0),
                           (y, ranges_log[0][1], ranges_log[0][2]), x_range,
                           annotations=[shift_arrow_base],
                           rendering_kw={"linestyles": "--", "linewidths": 1},
                           markers=[{'args': [new_cent, 0, 'x'], 'color':'red'}],
                           color='gray', **params_proj,
                           xlabel='y', ylabel='x')


# Собираем изображения в сетку общего рисунка
gs = GridSpec(2, 2)
mapping = {
    gs[0, 0]: p1 + p1_logos + p1_os, # + p1_log_legend,
    gs[0, 1]: p1_base_sh + p1_base,
    gs[1, 0]: p2 + p2_os + p2_logos,
    gs[1, 1]: p2_base_sh + p2_base,
}
plotgrid(gs=mapping, size=(7,7), legend=False);

Это можно сделать, например, задав переменное смещение для y-координаты.

Способ-1 (канонический)

Уравнение параболоида имеет вид $frac{x^2}{a^2} + frac{y^2}{b^2}+frac{z}{h}=0$ .

Такое мы использовали выше для задания формы зуба в разделе Зуб варана v.2 (параболоид).

Уравнение параболоида с смещением s по y имеет вид

$frac{x^2}{a^2} + frac{(y+s)^2}{b^2}+frac{z}{h}=0$

И если s — это просто некоторое положительное число, то весь параболоид сдвинется вдоль оси Oy влево. А если величина s зависит от z, то есть принимает какое-то свое значение для каждого горизонтального сечения параболоида, то мы имеем возможность сместить каждое сечение персонально на нужную нам величину. Схемы такого смещения показаны в правой части вышеприведенного рисунка.

Опираясь на знание, как ведет себя логарифм, и чувство прекрасного, зададим s как функцию от z следующего вида:

Тогда уравнение примет вид $frac{x^2}{a^2} + frac{(y+2ln(1-z))^2}{b^2}+frac{z}{h}=0$

Подставим значения a, b и h и получим искомое уравнение боковой поверхности искривленного драконьего зуба:

$frac{x^2}{4} + frac{(y+2ln(1-z))^2}{25}+frac{z}{25}=0$

Снизу, как обычно, объем зуба ограничен плоскостью .

Вводим полученное уравнение в сервисе desmos.com или аналогичном и наблюдаем изображение.

Способ 2 (параметрический)

Уравнение параболоида можно записать и в параметрическом виде, где аналогично способу-1 задать смещение s для y-координаты:

$begin{cases} z=-t^2 \ x=frac{a}{sqrt{h}} cdot t cdot cos(phi) \ y=frac{b}{sqrt{h}} cdot t cdot sin(phi) - s(z) end{cases}$

после подстановки найденной выше

$begin{cases} z=-t^2 \ x=frac{a}{sqrt{h}} cdot t cdot cos(phi) \ y=frac{b}{sqrt{h}} cdot t cdot sin(phi) - 2ln(1+t^2) end{cases}$

где параметры изменяются в следующих диапазонах: ,

Изобразим драконий зуб программно:

import numpy as np
import sympy as sp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from sympy.abc import x, y, z


a, b, h = 2, 5, 25

# функция нелинейного сдвига
log_line = 2*sp.log(1-z)

# векторизуем функцию логарифмического смещения для расчета сетки ниже
v_log_line = np.vectorize(lambda tz: float(log_line.subs({z:tz}).n()))

# сетка для переменных-параметров
tz = np.linspace(0, -np.sqrt(h), 30)
phi = np.linspace(0, 2 * np.pi, 30)  # одна сторона зуба

tz, phi = np.meshgrid(tz, phi)

# подставим сетку в параметрическое уравнение параболоида со смещением
Z = -tz**2
X = (a / np.sqrt(h)) * tz * np.cos(phi)
Y = (b / np.sqrt(h)) * tz * np.sin(phi) - v_log_line(-tz**2)

# График
ax = plt.figure().add_subplot(projection='3d')
ax.view_init(elev=30, azim=10, roll=0)

ax.plot_surface(X, Y, Z, color='w', edgecolor='royalblue', 
                rstride=2, cstride=2, alpha=0.3)

# Диапазоны и настройки графика
x_lims = [-2*a, 2*a]
y_lims = [-3*b, 0]
z_lims = [-h, 0]

ax.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
       xlabel='X', ylabel='Y', zlabel='Z',
       xticks=np.round(np.linspace(*x_lims, 5),1),
       yticks=np.round(np.linspace(*y_lims, 6),1)
       )
ax.set_aspect('equal')

plt.show()

Объем искривленного зуба варана

Разрежем мысленно фигуру на элементарные, предельно тонкие, слои — с высотой dh и площадью основания S. Объем каждого такого слоя будет равен :

Из параметрических уравнений поверхности . Найдем дифференциал

Площадь основания слоя S будет изменяться в зависимости от высоты, на которой вырезаем слой. По форме это в любом случае будет эллипс, но с различными полуосями. Обозначим их за a’ и b’ (штрихи просто как штрихи, не производные) и вспомним формулу площади эллипса:

Опять же из параметрических уравнений поверхности явно видны выражения для полуосей (множители перед синусом и косинусом):

$a'=frac{ acdot t}{sqrt{h}}$ , $b'=frac{ bcdot t}{sqrt{h}}$

Подставим их в формулу площади эллипса:

$S=pi cdot a' cdot b'=-pi frac{a cdot t}{sqrt{h}} frac{b cdot t}{sqrt{h}} 2t=2frac{pi a b}{h} t^3$

А ее, в свою очередь, подставим в выражение для элементарного объема

$dV=-2frac{pi a b}{h} t^3 dt$

Берем интеграл от левой и правой частей

$V=int_{-sqrt{h}}^0 frac{-2 pi a b}{h} t^3 dt$

Интегрируем мы в данной статье все численно, не будем делать исключение и для этого интеграла. Подставим значения a,b и h и найдем объем с помощью функции quad() библиотеки SciPy

from scipy import integrate
import numpy as np


a, b, h = 2, 5, 25

# Подынтегральная функция
def f(t):
 return - 2 * np.pi * a * b * t**3 / h

# Интеграл
V, err = integrate.quad(f, -np.sqrt(h), 0)

print(f'V = {V}')

V = 392.7 мм³.

Объем получился такой же, как для зуба в форме обычного параболоида. Что логично, поскольку мы лишь параллельно сдвинули слои в сторону, не изменив ни их форму, ни толщину. Аналогично тому, как если положить на стол колоду карт и сдвинуть верхние слои относительно нижних, объем колоды останется прежним (привет принципу Кавальери).

Объем искривленного зуба варана с полостью

В заключение учтем тот факт, что зуб варана имеет полость внутри [4]. Ориентируясь на продольный срез зуба, приведенный на Fig.3 из [4], опишем форму полости как аналогично параболическую с высотой 80% от высоты зуба и размерами основания в 30% от соответствующих размеров основания зуба. То есть с h = 20 мм, a = 0.6 мм, b = 1.5 мм. Подставив эти значения в формулу/код из предыдущего раздела получим V_{полости} = 28.3 мм³

И значит объем костной ткани зуба V = V_общий - V_{полости} = 392.7 - 28.3 = 364.4 мм³.

Кальций. Сколько вешать в граммах

Зная объем и плотность вещества можно найти его массу. Объем в наличии, осталось определится с веществом и плотностью.

Зуб варана состоит из дентина и покрывающей его эмали. Ввиду того, что слой эмали у комодского дракона чрезвычайно тонкий — 20 микрометров [3] - в рамках данной модели им пренебрегаем и считаем, что весь найденный выше объем занят дентином.

Каково экспериментальное значение плотности дентина у комодских варанов, мне, к сожалению, найти не удалось (если кто-то из читателей владеет такой информацией, буду признателен). Возьмем в качестве приближения параметры зубных тканей крокодила — рептилии схожей с вараном по размеру и хищному способу добывания пищи.

Минералом, придающим прочность зубам крокодила, является гидроксиапатит [5].

Найдем, каков процент его массы составляет чистый кальций.

Химическая формула гидроксиапатита [6], плотность =3.14 г/см³ [7] .

Массы атомов, входящих в формулу: $m_{Ca}=40.078 u$ , $m_{P}=30.974 u$ , $m_{O}=15.999 u$ , $m_{H}=1.008 u$ (где u — атомная единица массы).

Откуда общая масса рассматриваемой молекулы

$m_{total}=m_{Ca} cdot 5 + (m_{P} + m_{O} cdot 4) cdot 3 + (m_{O} + m_{H})=$

= 40.078 ∙ 5 + (30.974 + 15.999 ∙ 4) ∙ 3 + (15.999 + 1.008) = 502.307 u

И значит на долю кальция по массе приходится $eta_{mol}=(5 m_{Ca}) / m_{total}$ = (40.078 ∙ 5) / 502.307 = 0.3989 = 39.89 %

Содержание самого гидроксиапатита в дентине составляет $eta_{dent}$ = 0.713 = 71.3 % . Остальное — органические компоненты [5].

Таким образом, чистого кальция в дентине содержится $eta_{Ca}=eta_{dent} cdot eta_{mol}$ = 0.713 ∙ 0.3989 = 0.2844 = 28.44 % .

Оставшиеся 42.9% — это суммарная доля фосфора, кислорода и водорода.

Массовые доли рассматриваемых компонентов вещества зуба показаны на диаграмме ниже. Минеральная часть выделена сплошной заливкой, органическая — заполнением точками.

Найдем массу дентина одного драконьего зуба.

$m_{dent}=V cdot rho$ = 364.4 мм³ ∙ 3.14 ∙ 10^-3 г/мм³ = 1.144 г

Значит масса чистого кальция в одном зубе комодского варана составит

$m_{Ca}=eta_{Ca} cdot m_{dent}$ = 0.325 г

Поскольку всего зубов у дракона с острова Комодо 60 штук, для их полной замены необходимо 60 ∙ $m_{Ca}$ = 14.1 г кальция.

Если поделить это на 40 дней (период обновления зуба согласо [4]), то получим 0.353 гр./сутки = 353 мг/сутки.

Именно столько кальция требуется варану для обеспечения динамики регулярного обновления комплекса зубов.

Норма потребления по кальцию в целом будет выше, т.к. кроме зубов он содержится в костной ткани, необходим для сокращения скелетных мышц, передачи нервных импульсов и ряда других процессов в организме. К примеру, норма потребления кальция для человека в возрасте от 16 до 50 лет, когда роста зубов уже не происходит, составляет 800-1000 мг/сутки [8].

Режущая кромка

Режущая кромкаКромка v.1 (ровная)

Зуб комодского дракона замечателен не только своей формой, но и наличием укрепленной железом кромки [3]. На рисунке ниже она выделяется оранжевым цветом и волнистой формой.

Для начала опишем геометрию ровной кромки, без учета идущей по ней «волны».

Уравнение кривой и визуализация

Поскольку кромка зуба это линия, принадлежащая его поверхности, воспользуемся полученным ранее параметрическим представлением поверхности.

$begin{cases} z=-t^2 \ x=frac{a}{sqrt{h}} cdot t cdot cos(phi) \ y=frac{b}{sqrt{h}} cdot t cdot sin(phi) - 2ln(1+t^2) end{cases}$

Кромку будет образовывать множество точек со значениями углового параметра и (соответствующие противоположным краям), при том, что параметр t, отвечающий за положение по высоте, пробегает прежний диапазон значений .

Такое представление пригодится нам для трехмерной визуализации.

Однако в данном случае (и для будущего расчета длины кромки) можно обойтись и двумерным графиком в центральной вертикальной плоскости (), содержащей оба противоположных режущих края.

При фиксированном x = 0 в системе остается два уравнения:

$begin{cases} z=-t^2 \ y=frac{b}{sqrt{h}} cdot t cdot sin(phi) - 2ln(1+t^2) end{cases}$

Которые можно свести в одно, выразив t через z из первого уравнения и подставив во второе. Итого

$y=frac{b}{sqrt{h}} cdot sqrt{-z} cdot sin(phi)- 2ln(1-z)$

Подставим известные значения b и h.

Для переднего края (обозначен на рисунке ниже L₁) с $phi=frac{3pi}{2}$ приходим к следующему уравнению в декартовых координатах:

$y_1=-sqrt{-z} — 2ln(1-z)$

Для заднего края (обозначен на рисунке ниже L₂) с $phi=frac{3pi}{2}$ — к уравнению:

$y_2=sqrt{-z} — 2ln(1-z)$

Построим графики:

Код для построения графиков

import numpy as np
import sympy as sp
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from sympy.abc import x, y, z


a, b, h = 2, 5, 25

# Кривая смещения
log_line = 2*sp.log(1-z)

# Поле рисунка для построения графиков
fig = plt.figure(constrained_layout=True, figsize=(8,5))
spec = fig.add_gridspec(ncols=2, nrows=1)
spec_right = spec[1].subgridspec(ncols=1, nrows=5)

ax1 = fig.add_subplot(spec[0, 0], projection='3d')
ax2 = fig.add_subplot(spec_right[1:4])

# ------- График-1: зуб с выделенной режущей кромкой --------

# векторизуем функцию логарифмического смещения для расчета сетки ниже
v_log_line = np.vectorize(lambda tz: float(log_line.subs({z:tz}).n()))

# сетка для переменных-параметров
tz = np.linspace(0, np.sqrt(h), 30)
phi_0 = np.linspace(0, 2 * np.pi, 30)  # поверхность зуба
phi_1 = np.linspace((3/2)*np.pi, (3/2)*np.pi, 30)  # режущая кромка зуба
phi_2 = np.linspace(np.pi / 2, np.pi/ 2, 2)  # режущая кромка зуба

tz_0, phi_0 = np.meshgrid(tz, phi_0)
tz_1, phi_1 = np.meshgrid(tz, phi_1)
tz_2, phi_2 = np.meshgrid(tz, phi_2)

# График
ax1.view_init(elev=30, azim=-35, roll=0)

for tz, phi, col, lw in ([tz_0, phi_0, 'royalblue', 1], [tz_1, phi_1, 'r', 2], [tz_2, phi_2, 'r', 2]):
  # подставим сетку в параметрическое уравнение параболоида со смещением
  Z = -tz**2
  X = (a / np.sqrt(h)) * tz * np.cos(phi)
  Y = (b / np.sqrt(h)) * tz * np.sin(phi) - v_log_line(-tz**2)

  ax1.plot_surface(X, Y, Z, color='w', edgecolor=col, lw=lw,
                  rstride=5, cstride=5, alpha=0.3)

# Диапазоны и настройки графика
x_lims = [-2*a, 2*a]
y_lims = [-3*b, 0]
z_lims = [-h, 0]

ax1.set_aspect('equal')
ax1.set(xlim=x_lims, ylim=y_lims, zlim=z_lims,
        xlabel='X', ylabel='Y', zlabel='Z',
        xticks=np.round(np.linspace(*x_lims, 5),1),
        yticks=np.round(np.linspace(*y_lims, 6),1)
        )

# -------- График-2: режущая кромка на плоскости ---------

# Уравнение боковой поверхности изогнутого зуба
left_part_eq = x**2/a**2 + (y + log_line)**2/b**2 + z/h
equations = [left_part_eq]

sol1, sol2 = sp.solve(equations, y, dict=True)
print(sol1, sol2)

# Проекция вдоль Ox на плоскость Oyz
zy_section1 = sol1[y].subs({x:0})
zy_section2 = sol2[y].subs({x:0})

print(zy_section1, zy_section2)

# Диапазоны
zs = np.linspace(*z_lims, 30)
ys1 = [zy_section1.subs({z:val}).n() for val in zs]
ys2 = [zy_section2.subs({z:val}).n() for val in zs]

# График
ax2.plot(ys1, zs, label='$L_1$')
ax2.plot(ys2, zs, label='$L_2$')

# Настройки графика
ax2.set(xlim=y_lims, ylim=z_lims, xlabel='Y', ylabel='Z')
ax2.set_aspect('equal')
ax2.grid()
ax2.legend()

plt.show()

Длина ровной режущей кромки

Длину кромки найдем при помощи криволинейного интеграла 1го рода.

Длина передней кромки:

$l_1=intlimits_{L_1} dl=intlimits_{-h}^0 sqrt{1+frac{dy_1}{dz}} dz=intlimits_{-h}^0 sqrt{1-frac{2}{1 — z} - frac{sqrt{-z}}{2z}} dz$

Длина задней кромки вычисляется по той же схеме, отличаясь лишь знаками в выражении для :

$l_2=intlimits_{-h}^0 sqrt{1+frac{dy_2}{dz}} dz=intlimits_{-h}^0 sqrt{1-frac{2}{1 - z} + frac{sqrt{-z}}{2z}} dz$

И полная длина

Возьмем данные интегралы численно и определим итоговую длину.

import sympy as sp 
from sympy.abc import x, y, z 
from scipy import integrate 


a, b, h = 2, 5, 25

# Кривая смещения оси 
log_line = 2*sp.log(1-z)

# Уравнение боковой поверхности изогнутого зуба 
left_part_eq = x**2/a**2 + (y + log_line)**2/b**2 + z/h 
equations = [left_part_eq] 
# разрешим относительно y
sol1, sol2 = sp.solve(equations, y, dict=True) 
print('Решения отностительно y: n', sol1, sol2)

# Проекция вдоль Ox на плоскость Oyz 
zy_section1 = sol1[y].subs({x:0}) 
zy_section2 = sol2[y].subs({x:0}) 
print('Проекции: n', zy_section1, zy_section2)

# Берем производные 
derivative_y1 = sp.diff(zy_section1, z) 
derivative_y2 = sp.diff(zy_section2, z) 
print('Их производные: n', derivative_y1, derivative_y2)


# Подынтегральная функция для криволинейного интеграла 
def f(z_val, curve_name):
  derivative = {'L1': derivative_y1, 'L2': derivative_y2}

  return (1 + derivative[curve_name].subs({z:z_val}).n()**2) ** 0.5


# Интеграл 
length = integrate.quad(f, -h, 0, 'L1')[0] + integrate.quad(f, -h, 0, 'L2')[0] 
print(f'Длина кромки {length} мм')

Таким образом длина ровного режущего края l=54.18 мм

Кромка v.2 (волнистая)

Циклоида путешественница

Для начала вспомним, как выглядит параметрическое уравнение циклоиды:

$begin{cases} x=R(t-sin(t)) \ y=R(1-cos(t)) end{cases}$

Если производящая окружность будет катиться не по горизонтальной прямой , а по кривой , то уравнение примет следующий вид:

$begin{cases} x=R(t-sin(t)) \ y=R(1-cos(t))+f(x(t)) end{cases}$

Как можно заметить, здесь использован тот же принцип с заданием функции смещения по y, что и при искривлении оси параболоида в разделе Зуб варана v.3 (логарифм).

Такая кривая является частным случаем квазитрохоиды (которая в общем случае описывает сочетание произвольного числа разновидностей поступательного и вращательного движений).

Демонстрационный пример

Рассмотрим демонстрационный пример — с опорной кривой вида , по которой отправим в путешествие производящую окружность радиусом R=0.5.

$begin{cases} x=0.5(t-sin(t)) \ y=0.5(1-cos(t))+0.1x(t) end{cases}$

или

$begin{cases} x=0.5(t-sin(t)) \ y=0.5(1-cos(t)) + 0.05 (t-sin(t)) end{cases}$

Что на графике будет выглядеть так:

Код для построения графика:

import matplotlib.pyplot as plt 
import numpy as np 


R = 0.5  # радиус производящей окружности

# Уравнение циклоиды, опирающейся на наклонную прямую 
x = lambda t: R * (t — np.sin(t))
line = lambda t: 0.1 * x(t)  # опорная прямая 
y = lambda t: R * (1 — np.cos(t)) + line(t)

# Подготовка данных для графика 
ts = np.linspace(0, 6*np.pi, 100) 
xs = [x(t) for t in ts] 
ys = [y(t) for t in ts] 
framework = [line(t) for t in ts]

# График 
fig, ax = plt.subplots()
ax.plot(xs, ys, label='L') 
ax.plot(xs, framework)

# Оформление графика
ax.grid() 
ax.legend()

plt.show()

Уравнение волнистой режущей кромки и визуализация

Теперь направим производящую окружность квазитрохоиды вдоль режущей кромки зуба, уравнения для которой мы получили ранее в разделе Кромка v.1 (ровная).

Передняя кромка зуба описывается уравнением $y_1=-sqrt{-z} - 2ln(1-z)$

Значит параметрическое уравнение квазитрохоиды, опирающейся на эту кромку будет иметь вид

$begin{cases} z=0.5(t-sin(t)) \ y=0.5(1-cos(t))+y_1(z(t)) end{cases}$

Подставим выражение для y₁

$begin{cases} z=0.5(t-sin(t)) \y=0.5(1-cos(t)) - sqrt{-z(t)} - 2ln(1-z(t)) end{cases}$

Заменим z во втором уравнении на ее выражение через параметр t из первого уравнения

$begin{cases} z=0.5(t-sin(t)) \ y=0.5(1-cos(t)) - sqrt{0.5(sin(t)-t)} - 2ln(1-0.5(t-sin(t))) end{cases}$

Получили параметрическое уравнение волнистой передней кромки зуба при радиусе производящей окружности R=0.5.

Такое значение радиуса обеспечит нам наглядность чертежа. Далее при расчете длины возьмем значение по экспериментальным данным (меньше на порядок).

Аналогично и для задней кромки. Уравнение опорной кривой имеет вид $y_2=sqrt{-z} - 2ln(1-z)$

И, значит, параметрическое уравнение для задней волнистой кромки:

$begin{cases} z=0.5(t-sin(t)) \ y=0.5(1-cos(t)) + sqrt{0.5(sin(t)-t)} - 2ln(1 - 0.5t + 0.5 sin(t))) end{cases}$

Параметр t, пробегает значения от до 0, где — решение уравнения .

Воспроизведем это в программном коде и графической форме.

Первым этапом находим уравнения для ровных передней и задней кромок.

import matplotlib.pyplot as plt 
import numpy as np 
import sympy as sp 
from sympy.abc import x, y, z, u 

a, b, h = 2, 5, 25
R = 0.5  # радиус производящей окружности

# Кривая смещения оси
log_line = 2*sp.log(1-z)

# Уравнение боковой поверхности изогнутого зуба
left_part_eq = x2/a2 + (y + log_line)2/b2 + z/h
equations = [left_part_eq]
# Разрешим относительно y
sol1, sol2 = sp.solve(equations, y, dict=True)
print('Решения относительно y: n', sol1, sol2)

# Проекция вдоль Ox на плоскость Oyz
zy_section1 = sol1[y].subs({x:0})
zy_section2 = sol2[y].subs({x:0})
print('Проекции: n', zy_section1, zy_section2)

Получили выражения для передней и задней ровных кромок (в коде выше zy_section1 и zy_section2). Применим их в качестве направляющих для развертывания волнообразного режущего края.

# Параметрическое уравнение квазитрохоиды на направляющих кромках
Z = lambda t: R * (t - np.sin(t))

# Направляющие
line_1 = lambda t: float(zy_section1.subs({z:Z(t)}).n())
line_2 = lambda t: float(zy_section2.subs({z:Z(t)}).n())

Y_1 = lambda t: -R * (1 - np.cos(t)) + line_1(t)
Y_2 = lambda t: R * (1 - np.cos(t)) + line_2(t)

Определим пределы изменения переменных

# Границы и диапазоны
y_lims = [-3*b, 1]
z_lims = [-h, 1]

# при каком значении параметра z = -h (т.е. на нижнем крае)
equation = sp.Eq(R * (u - sp.sin(u)), -h)
bottom_u = float(sp.nsolve(equation, u, -1))

# z = 0 очевидно при u = 0
top_u = 0

print(f'Пределы изменения параметра u: от {bottom_u} до {top_u}')

# генерируем данные для построения графика
ts = np.linspace(bottom_u, top_u, 100)

ys_1, ys_2 = [Y_1(t) for t in ts], [Y_2(t) for t in ts]
zs_1, zs_2 = [Z(t) for t in ts], [Z(t) for t in ts]
framework_y1 = [float(zy_section1.subs({z:z_val}).n()) for z_val in zs_1]
framework_y2 = [float(zy_section2.subs({z:z_val}).n()) for z_val in zs_2]

Построим график

# График
fig, ax = plt.subplots()

ax.plot(ys_1, zs_1, label='$C_1$')
ax.plot(framework_y1, zs_1, '--', label='$L_1$')
ax.plot(ys_2, zs_2, label='$C_2$')
ax.plot(framework_y2, zs_2, '--', label='$L_2$')

# параметры отображения графика
ax.set(xlim=y_lims, ylim=z_lims, xlabel='Y', ylabel='Z')
ax.set_aspect('equal')
ax.grid()
ax.legend()

plt.show()

Длина волнистой режущей кромки

Длина кривой L, заданной параметрически в зависимостями и , где , будет вычисляться по формуле

$l=intlimits_{L}dl=intlimits_{t_1}^{t_2} sqrt{(z'_t)^2+(y'_t)^2} dt$

Займемся ее вычислением.

Зададимся уравнением боковой поверхности зуба и найдем уравнения для передней и задней кромок. Эту часть кода мы уже использовали в предыдущих разделах. Единственным отличием здесь будет значение радиуса порождающей окружности. Установим его равным R=0.01 мм в соответствии с экспериментальными данными [3]

Уравнения волнистых режущих кромок (Sympy)

import sympy as sp 
from sympy.abc import x, y, z, u 
from scipy import integrate 

a, b, h, R = 2, 5, 25, 0.01

# Кривая смещения оси 
log_line = 2*sp.log(1-z)

# Уравнение боковой поверхности изогнутого зуба 
left_part_eq = x**2/a**2 + (y + log_line)**2/b**2 + z/h 
equations = [left_part_eq] 
# Разрешим относительно y
sol1, sol2 = sp.solve(equations, y, dict=True) 
print('Решения отностительно y: n', sol1, sol2)

# Проекция вдоль Ox на плоскость Oyz 
zy_section1 = sol1[y].subs({x:0}) 
zy_section2 = sol2[y].subs({x:0}) 
print('Проекции: n', zy_section1, zy_section2

Используем полученные как проекции выражения для передней и задней ровных кромок в качестве направляющих для развертывания волнообразного режущего края

# Параметрические уравнения передней и задней кромки 
zu = R * (u — sp.sin(u)) 

# Опорные кривые
line_1 = zy_section1.subs({z:zu}) 
line_2 = zy_section2.subs({z:zu}) 

Y_1 = -R * (1 — sp.cos(u)) + line_1 
Y_2 = R * (1 — sp.cos(u)) + line_2

После чего посчитаем необходимые производные и возьмем интеграл

# Их производные по параметру 
derivative_zu = sp.diff(zu, u) 
derivative_yu1 = sp.diff(Y_1, u) 
derivative_yu2 = sp.diff(Y_2, u) 
print('Их производные: n', derivative_zu, derivative_yu1, derivative_yu2) 


# Подынтегральная функция для криволинейного интеграла 
def f(u_val, curve_name):
  derivative = {'L1': derivative_yu1, 'L2': derivative_yu2}

  return (derivative_zu.subs({u:u_val}).n()**2 +
          derivative[curve_name].subs({u:u_val}).n()**2) ** 0.5

  
# Пределы 
# при каком значении параметра z = -h (т.е. на нижнем крае) 
equation = sp.Eq(R * (u — sp.sin(u)), -h) 
bottom_u = float(sp.nsolve(equation, u, -1))

# z = 0 очевидно при u = 0 
top_u = 0

# Интеграл
length = integrate.quad(f, bottom_u, top_u, 'L1')[0] + integrate.quad(f, bottom_u, top_u, 'L2')[0]

print(f'Длина кромки {length} мм')

Итого, длина волнообразной режущей кромки составляет 66.47 мм. Циклоидные дуги увеличили длину кромки на 23%.

Объем волнистой режущей кромки

Из [3] можно заключить, что характерный размер идущего по кромке оранжевого канта в поперечном сечении составляет 20 мкм. Представим его как цилиндр диаметром D = 20 мкм, уложенный по найденной выше линии волнистой кромки и рассчитаем объем.

Радиус сечения r = D / 2 = 10 мкм = 0.01 мм.

Площадь сечения

Объем цилиндра

Итак, объем железосодержащего "шнура", образующего режущую кромку, равен V = 0.021 мм³.

Железо оранжевой кромки

Согласно исследованию [3] оранжевая режущая кромка зуба варана состоит из ферригедрита. Формула ферригидрита: [9] Плотность =3.96 г/см³ = 3.96∙10^-3 г/мм³ [10]

Массы атомов, входящих в формулу: $m_{Fe} =55.845 u$ , , (где u — атомная единица массы).

Откуда общая масса рассматриваемой молекулы $m_{total}=5 (2 m_{Fe} + 3 m_O) + 9 (2 m_H + m_O)$ = 5 (55.845 ∙ 2 + 15.999 ∙ 3) + 9 * (1.008 ∙ 2 + 15.999) = 960.57 u

Массовая доля Fe составляет

$eta_{Fe}=(10 m_{Fe}) / m_{total}$ = (10 ∙ 55.845) / 960.57 = 0.581 = 58.1 %

Масса чистого железа в покрытии режущей кромки зуба комодского варана составит

$m_{Fe}=eta_{Fe} cdot rho cdot V$ = 0.581 ∙ 3.96∙10^-3 ∙ 0.021 = 0.048 ∙ 10^-3 г = 0.048 мг

Масса железа на кромках полного комплекта из 60 зубов $m_{teeth}=60 cdot m_{Fe}$ = 60 ∙ 0.048 мг = 2.88 мг

Зная, что период замены зуба у варана составляет 40 суток [4], найдем, сколько железа нужно в день при условии его равномерного расходования на покрытие формирующихся режущих кромок. $m_{teeth}$ / 40 = 2.88 мг/ 40 суток = 0.072 мг/сутки.

Не так уж много. Если сравнить с человеком, то средняя суточная доза потребления железа для взрослых мужчин и пожилых людей обоих полов составляет 8 мг (для женщин в 2-3 раза больше) [11]. При этом оно расходуется в основном на не связанные напрямую с зубами нужды организма.

Варану на покрытие острых кромок своих зубов нужно порядка 1% от такого объема.

Заключение

В рамках рассмотренных моделей мы получили, что варану для обеспечения динамики обновления зубов необходимо порядка 353 мг кальция и 72 мкг железа в сутки. Для сравнения: это примерно половина суточной потребности человека в кальции и менее одного процента суточной потребности человека в железе соответственно. Цифры не астрономические. И для комодского дракона, с его охотничьими способностями и положением на вершине пищевой цепи своего региона, вполне достижимые.

Что касается зубов человека, биоинженерные технологии понемногу идут в сторону создания (включения) механизма выращивания новых зубов взамен поврежденных или отсутствующих ([12] или подробнее в [13]). Но, думаю, еще долгое время варану конкуренция в этом вопросе с нашей стороны не грозит.

Источники

Автор: shbma

Источник

Информация

Обсуждаемое

Рекомендуем

Зуб комодского дракона

Оглавление

Введение

Геометрия зуба

Зуб варана v.1 (конус)

Объем конического зуба

Уравнение поверхности и визуализация

Зуб варана v.2 (параболоид)

Уравнение поверхности и визуализация

Объем параболического зуба

Зуб варана v.3 (логарифм)

Уравнение поверхности и визуализация

Объем искривленного зуба варана

Объем искривленного зуба варана с полостью

Кальций. Сколько вешать в граммах

Режущая кромка

Режущая кромкаКромка v.1 (ровная)

Уравнение кривой и визуализация

Длина ровной режущей кромки

Кромка v.2 (волнистая)

Циклоида путешественница

Уравнение волнистой режущей кромки и визуализация

Длина волнистой режущей кромки

Объем волнистой режущей кромки

Железо оранжевой кромки

Заключение

Источники

Архив

Информация

Обсуждаемое

Рекомендуем

Зуб комодского дракона

Оглавление

Введение

Геометрия зуба

Зуб варана v.1 (конус)

Объем конического зуба

Уравнение поверхности и визуализация

Зуб варана v.2 (параболоид)

Уравнение поверхности и визуализация

Объем параболического зуба

Зуб варана v.3 (логарифм)

Уравнение поверхности и визуализация

Объем искривленного зуба варана

Объем искривленного зуба варана с полостью

Кальций. Сколько вешать в граммах

Режущая кромка

Режущая кромкаКромка v.1 (ровная)

Уравнение кривой и визуализация

Длина ровной режущей кромки

Кромка v.2 (волнистая)

Циклоида путешественница

Уравнение волнистой режущей кромки и визуализация

Длина волнистой режущей кромки

Объем волнистой режущей кромки

Железо оранжевой кромки

Заключение

Источники

Рекомендованный контент

Новости

Актуальные темы

Архив